Re: [petsc-users] MatPtAP

[Samuel Lanthaler]

Ok, I now realize that I had implemented the boundary conditions in an 
unnecessarily complicated way... As you pointed out, I can just 
manipulate individual matrix rows to enforce the BC's. In that way, I 
never have to call MatPtAP, or do any expensive operations. Probably 
that's what is commonly done. I've changed my code, and it seems to work 
fine and is much faster.


Thanks a lot for your help, everyone! This has been very educational for me.

Samuel


On 06/01/2018 09:04 PM, Hong wrote:

Samuel,
I have following questions:
1) Why solving \lambda*Pt*B*P*y=Pt*A*P*y is better than solving 
original \lambda*B*x = A*x?


2) Does your eigen solver require matrix factorization? If not, i.e., 
only uses mat-vec multiplication, then you may implement
z = (Pt*(A*(P*y))) in your eigensolver instead of using mat-mat-mat 
multiplication.


3) petsc PtAP() was implemented for multigrid applications, in which 
the product C = PtAP is a denser but a much smaller matrix.
I have not seen the use of your case, that P is square with same size 
as A. If C is much denser than A, then PtAP consumes a large portion 
of time is anticipated.


Hong

On Fri, Jun 1, 2018 at 12:35 PM, Smith, Barry F. <mailto:bsm...@mcs.anl.gov>> wrote:



  Could you save the P matrix with MatView() using a binary viewer
and the A matrix with MatView() and the binary viewer and email
them to petsc-ma...@mcs.anl.gov <mailto:petsc-ma...@mcs.anl.gov> ?
Then we can run the code in the profiler with your matrices and
see if there is any way to speed up the computation.

   Barry


> On Jun 1, 2018, at 11:07 AM, Samuel Lanthaler
mailto:s.lantha...@gmail.com>> wrote:
>
> On 06/01/2018 03:42 PM, Matthew Knepley wrote:
>> On Fri, Jun 1, 2018 at 9:21 AM, Samuel Lanthaler
mailto:s.lantha...@gmail.com>> wrote:
>> Hi,
>>
>> I was wondering what the most efficient way to use MatPtAP
would be in the following situation: I am discretizing a PDE
system. The discretization yields a matrix A that has a band
structure (with k upper and lower bands, say). In order to
implement the boundary conditions, I use a transformation matrix P
which is essentially the unit matrix, except for the entries
P_{ij} where i,j>
>> P =  [ B, 0, 0, 0, ..., 0, 0 ]
>>[  0, 1, 0, 0, ..., 0, 0 ]
>>[  ]
>>[  ]
>>[  ..., 1, 0 ]
>>[  0, 0, 0, 0, ..., 0, C ]
>>
>> with B,C are (k-by-k) matrices.
>> Right now, I'm simply constructing A, P and calling
>>
>> CALL

MatPtAP(petsc_matA,petsc_matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,petsc_matPtAP,ierr)
>>
>> where I haven't done anything to pestc_matPtAP, prior to this
call. Is this the way to do it?
>>
>> I'm asking because, currently, setting up the matrices A and P
takes very little time, whereas the operation MatPtAP is taking
quite long, which seems very odd... The matrices are of type
MPIAIJ. In my problem, the total matrix dimension is around 10'000
and the matrix blocks (B,C) are of size ~100.
>>
>> Are you sure this is what you want to do? Usually BC are local,
since by definition PDE are local, and
>> are applied pointwise. What kind of BC do you have here?
>>
>
> The boundary conditions are a mixture of Dirichlet and Neumann;
in my case, the PDE is a system involving 8 variables on a disk,
where the periodic direction is discretized using a Fourier series
expansion, the radial direction uses B-splines.
>
> In reality, I have two matrices A,B, and want to solve the
eigenvalue problem \lambda*B*x = A*x.
> I found it quite convenient to use a transformation P to a
different set of variables y, such that x=P*y and x satisfies the
BC iff certain components of y are 0. The latter is enforced by
inserting spurious eigenvalues at the relevant components of y in
the transformed eigenvalue problem \lambda*Pt*B*P*y=Pt*A*P*y.
After solving the EVP in terms of y, I get back x=P*y.
> Is this an inherently bad/inefficient way of enforcing BC's? Thanks.
>
>
>
>
>>   Thanks,
>>
>> Matt
>>
>> Thanks in advance for any ideas.
>>
>> Cheers,
>> Samuel
>>
>>
>>
>> --
>> What most experimenters take for granted before they begin
their experiments is infinitely more interesting than any results
to which their experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
<https://www.cse.buffalo.edu/%7Eknepley/>
>

Re: [petsc-users] MatPtAP

[Samuel Lanthaler]

On 06/01/2018 04:55 PM, Smith, Barry F. wrote:

How do you know it is exactly the MatPtAP() routine that is taking the 
large time?


At least according to a simple timing (CALL CPU_TIME(...)), the call to 
MatPtAP appears to consume more time than the actual discretization of 
the PDE. Let me attach the output of -log_view.




Please send the output of a run with -log_view so we can see where the time 
is being spent.

Barry



On Jun 1, 2018, at 8:21 AM, Samuel Lanthaler  wrote:

Hi,

I was wondering what the most efficient way to use MatPtAP would be in the following 
situation: I am discretizing a PDE system. The discretization yields a matrix A that 
has a band structure (with k upper and lower bands, say). In order to implement the 
boundary conditions, I use a transformation matrix P which is essentially the unit 
matrix, except for the entries P_{ij} where i,j

 *  VENUS-LEVIS: rev Unversioned directory *
 --
  DELTA-F COMPUTATION FOR 
HYBRID KINETIC-MHD
 --
 **
 Cleaning sub-directories and creating diagnostic tree
 ---
 . reading simulation data
 . reading equilibrium data
 . reading IO specification
 synchronize nparts,nprocs!
 nparts must be multiple of nprocs
 nparts,nprocs   1   4
 Resetting values to satisfy constraint:
 nparts,nprocs   4   4
 . initializing magnetic equilibrium
 Loading SFL (STRAIGHT-FIELD LINE) equilibrium
 - Plasma volume   197.391538791518 
  -- 
 | magnetic equilibrium: 
 | - 
 | B0 =   0.984245824471487 
 | R0 =10.1237829388870 
 | alfven velocity  =4809804.25531587 
 | alfven frequency =475099.504241706 
  -- 
 
 Profiles related to flow (on axis):
 R0: 10.1237829388870 
 mOmp:   2.00 
 T0 [keV]:   1.1099122785 
 Omega0 [Hz]:   0.000E+000
 Mach0: 0.000E+000
 
 . initializing analytic distributions
 . bcasting analytic profiles data
 . calling test_analytic_profiles
 . reading equilibrium distribution data
 . bcasting equilibrium distribution data
 - TAE updaters
 . dumping analytic distribution; Ne,Np,Nmu  100 100
   5
 - TAE updaters
 . initializing delta-f computation
 . reading delta-f grid data
 . bcasting delta-f grid data
 -
 WARNING: Using df_coll_freq to set ftime%dE_prefix.
 ftime%dE_prefix (original):  (0.000E+000,0.000E+000)
 ftime%dE_prefix (set):   (4.75099504241706,237.549752120853)
 -
 ftime%grate =4.75099504241706 
 ftime%freq =237.549752120853 
 ftime%dE_prefix =  (4.75099504241706,237.549752120853)
 -
 ftoro%n  =1
 ftoro%dE_prefix =  (0.000E+000,-1.00)
 *   STARTING SIMULATION   *
 over4  processors
 with4  orbits
 **
 initializing from MHD in
 ... setting type
 PDE.f90--> Chosen case: 
 *   MHDflow
 ... initializing elements
 initializing with bunching on rational surfaces:
 q =1.00 
 q =2.00 
 --
 
 q = qvals found for:
 q =1.1265158867 at s =   0.27216000131 
 q =2.0720104897 at s =   0.7999601 
 
 Nrad =  180
 mmin =   -3
 mmax =5
 ntor =   -1
 test_mode =7
 Initializing operator
 setting natural units in MHDflow8
 setting natural units in MHDflow8
 setting natural units in MHDflow8
 setting natural units in MHDflow8
 Discretizing operator
 Discretization in progress (computation of matrices)
 Assembly of matA,matB done.
 Accounting for use of mixed finite elements.
 Add boundary conditions
 Adding boundary conditions: Neumann/Dirichlet
 ... setting up BC
 ... Setting up R/T (transformation) matrices
 ... Applying BC to A/B
 ... A->RtAR, B->RtBR
 ---
 time MatPtAP:   17.39000 
 time MatZeroRowsColumns:   0.670 
 time MatAssembly:  4.0001336E-003
 -
 ---
 time MatPtAP:   16.35700 
 time MatZeroRowsColumns:   0.6319998 
 time MatAssembly:  3.114E-003
 -
 ... writing out matrices in boundary conditions...
 destroying local objects
 time set BC:43.90700 
 Discretization done.
 writing out matrices
 Writing 

Re: [petsc-users] MatPtAP

[Samuel Lanthaler]

On 06/01/2018 03:42 PM, Matthew Knepley wrote:
On Fri, Jun 1, 2018 at 9:21 AM, Samuel Lanthaler 
mailto:s.lantha...@gmail.com>> wrote:


Hi,

I was wondering what the most efficient way to use MatPtAP would
be in the following situation: I am discretizing a PDE system. The
discretization yields a matrix A that has a band structure (with k
upper and lower bands, say). In order to implement the boundary
conditions, I use a transformation matrix P which is essentially
the unit matrix, except for the entries P_{ij} where i,jAre you sure this is what you want to do? Usually BC are local, since 
by definition PDE are local, and

are applied pointwise. What kind of BC do you have here?



The boundary conditions are a mixture of Dirichlet and Neumann; in my 
case, the PDE is a system involving 8 variables on a disk, where the 
periodic direction is discretized using a Fourier series expansion, the 
radial direction uses B-splines.


In reality, I have two matrices A,B, and want to solve the eigenvalue 
problem \lambda*B*x = A*x.
I found it quite convenient to use a transformation P to a different set 
of variables y, such that x=P*y and x satisfies the BC iff certain 
components of y are 0. The latter is enforced by inserting spurious 
eigenvalues at the relevant components of y in the transformed 
eigenvalue problem \lambda*Pt*B*P*y=Pt*A*P*y. After solving the EVP in 
terms of y, I get back x=P*y.

Is this an inherently bad/inefficient way of enforcing BC's? Thanks.





  Thanks,

Matt

Thanks in advance for any ideas.

Cheers,
Samuel




--
What most experimenters take for granted before they begin their 
experiments is infinitely more interesting than any results to which 
their experiments lead.

-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/%7Emk51/>

[petsc-users] MatPtAP

[Samuel Lanthaler]

Hi,

I was wondering what the most efficient way to use MatPtAP would be in 
the following situation: I am discretizing a PDE system. The 
discretization yields a matrix A that has a band structure (with k upper 
and lower bands, say). In order to implement the boundary conditions, I 
use a transformation matrix P which is essentially the unit matrix, 
except for the entries P_{ij} where i,j

P =  [ B, 0, 0, 0, ..., 0, 0 ]
   [  0, 1, 0, 0, ..., 0, 0 ]
   [  ]
   [  ]
   [  ..., 1, 0 ]
   [  0, 0, 0, 0, ..., 0, C ]

with B,C are (k-by-k) matrices.
Right now, I'm simply constructing A, P and calling

CALL 
MatPtAP(petsc_matA,petsc_matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,petsc_matPtAP,ierr)


where I haven't done anything to pestc_matPtAP, prior to this call. Is 
this the way to do it?


I'm asking because, currently, setting up the matrices A and P takes 
very little time, whereas the operation MatPtAP is taking quite long, 
which seems very odd... The matrices are of type MPIAIJ. In my problem, 
the total matrix dimension is around 10'000 and the matrix blocks (B,C) 
are of size ~100.


Thanks in advance for any ideas.

Cheers,
Samuel

Re: [petsc-users] configuration option for PETSc on cluster

[Samuel Lanthaler]

Dear Satish,
Your workaround worked! Great!
Thanks a lot for all your help, everyone. =)

Cheers,
Samuel

On 04/24/2018 06:33 PM, Satish Balay wrote:

Hm - I'm not sure why configure didn't try and set -fPIC compiler options.

Workarround is:

CFLAGS=-fPIC FFLAGS=-fPICC CXXFLAGS=-fPIC

or alternatively - you can disable sharedlibraries

--with-shared-libraries=0

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Dear Satish,

As you suggested, the configure appears to work with

./configure --with-fc=ftn --with-cc=cc --with-cxx=CC --with-debugging=0
--with-scalar-type=complex --download-scalapack --download-mumps=yes
--download-superlu --with-batch --known-mpi-shared-libraries=1
--with-blaslapack-dir=${MKLROOT} --known-64-bit-blas-indices=0
--LDFLAGS="-dynamic"

However, after doing the reconfigure, I tried to compile using make, and ran
into trouble with CLINKER:

 [...]
 CC arch-complex/obj/sys/classes/draw/impls/x/xops.o
  CLINKER arch-complex/lib/libpetsc.so.3.09.0
Warning:
-dynamic was already seen on command line, overriding with -shared.
/usr/bin/ld: arch-complex/obj/sys/fsrc/somefort.o: relocation R_X86_64_32
against `petscfortran5_' can not be used when making a shared object;
recompile with -fPIC
arch-complex/obj/sys/fsrc/somefort.o: error adding symbols: Bad value
[...]
**ERROR*
   Error during compile, check arch-complex/lib/petsc/conf/make.log
   Send it and arch-complex/lib/petsc/conf/configure.log to
petsc-ma...@mcs.anl.gov

makefile:30: recipe for target 'all' failed
make: *** [all] Error 1

Any ideas? Again, I'm attaching the .log files.



On 04/24/2018 04:14 PM, Satish Balay wrote:

Sorry I should have clarified - its a configure option.

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Sorry, I'm confused: Can I add this somehow when doing the ./configure? Or
do
you mean to add it when compiling ex5f with the configuration I had
initially?

Cheers,
Samuel


On 04/24/2018 04:03 PM, Satish Balay wrote:

Ok - I do see: libmkl_intel_lp64.so libmkl_sequential.so libmkl_core.so -
these should
be the ones that get picked up.

I do not know why 'cc' is picking up static libraries instead of the
shared ones. Perhaps it needs an extra option

'LDFLAGS=-dynamic'

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Dear Satish,

I get the following when doing the ls:

[lanthale@daint103 mkl]$ ls
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/
libmkl_ao_worker.so libmkl_blacs_sgimpt_lp64.so
libmkl_intel_lp64.alibmkl_sequential.so
libmkl_avx2.so  libmkl_blas95_ilp64.a
libmkl_intel_lp64.so
libmkl_tbb_thread.a
libmkl_avx512_mic.solibmkl_blas95_lp64.a
libmkl_intel_thread.a
libmkl_tbb_thread.so
libmkl_avx512.solibmkl_cdft_core.a libmkl_intel_thread.so
libmkl_vml_avx2.so
libmkl_avx.so   libmkl_cdft_core.so
libmkl_lapack95_ilp64.a
libmkl_vml_avx512_mic.so
libmkl_blacs_intelmpi_ilp64.a   libmkl_core.a libmkl_lapack95_lp64.a
libmkl_vml_avx512.so
libmkl_blacs_intelmpi_ilp64.so  libmkl_core.so libmkl_mc3.so
libmkl_vml_avx.so
libmkl_blacs_intelmpi_lp64.alibmkl_def.so libmkl_mc.so
libmkl_vml_cmpt.so
libmkl_blacs_intelmpi_lp64.so   libmkl_gf_ilp64.a libmkl_pgi_thread.a
libmkl_vml_def.so
libmkl_blacs_openmpi_ilp64.alibmkl_gf_ilp64.so libmkl_pgi_thread.so
libmkl_vml_mc2.so
libmkl_blacs_openmpi_ilp64.so   libmkl_gf_lp64.a libmkl_rt.so
libmkl_vml_mc3.so
libmkl_blacs_openmpi_lp64.a libmkl_gf_lp64.so
libmkl_scalapack_ilp64.a
libmkl_vml_mc.so
libmkl_blacs_openmpi_lp64.solibmkl_gnu_thread.a
libmkl_scalapack_ilp64.so
locale
libmkl_blacs_sgimpt_ilp64.a libmkl_gnu_thread.so
libmkl_scalapack_lp64.a
libmkl_blacs_sgimpt_ilp64.solibmkl_intel_ilp64.a
libmkl_scalapack_lp64.so
libmkl_blacs_sgimpt_lp64.a  libmkl_intel_ilp64.so libmkl_sequential.a


Cheers,
Samuel

On 04/24/2018 03:50 PM, Satish Balay wrote:

Executing: cc  -o /tmp/petsc-tyd_Hm/config.libraries/conftest -O
/tmp/petsc-tyd_Hm/config.libraries/conftest.o
-Wl,-rpath,/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64
-L/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64
-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lstdc++ -ldl
Possible ERROR while running linker: exit code 256
stderr:
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(mkl_semaphore.o):
In function `mkl_serv_load_inspector':
mkl_semaphore.c:(.text+0x123): warning: Using 'dlopen' in statically
linked
applications requires at runtime the shared libraries from the glibc
version
used for linking
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(dgetrs.o):
In function `mkl_lapack_dgetrs':
dgetrs_gen.f:(.text+0x224): undefined reference to `mkl_blas_dtrsm'
dgetrs_gen.f:(.text+0x2b1): undefined

Re: [petsc-users] configuration option for PETSc on cluster

[Samuel Lanthaler]

... wait a minute, I found this on the cluster website

   Regardless of which programming environment (compiler suite) you are
   using, you should always compile code with the Cray wrapper commands
   (|cc| for C code, |CC| for C++ code and |ftn| for Fortran code).
   _The compiler wrappers produce statically-linked executables by
   default; for dynamic libraries use the __|-dynamic|__flag or
   __|export CRAYPE_LINK_TYPE=dynamic|__before building._ Note that
   dynamic linking is the default when the cudatoolkit module is
   loaded. You should include the compiler flags appropriate to the
   underlying compiler suite (Intel, GNU, PGI, Cray). There are also
   specific options for the
   wrappers themselves: see the man pages of |ftn|, |cc|, and |CC| for
   details

|| That sounds like what you were saying Satish. I'll try again with 
"|export CRAYPE_LINK_TYPE=dynamic"| and come back to you! Thanks.




On 04/24/2018 04:03 PM, Satish Balay wrote:

Ok - I do see: libmkl_intel_lp64.so libmkl_sequential.so libmkl_core.so - these 
should
be the ones that get picked up.

I do not know why 'cc' is picking up static libraries instead of the
shared ones. Perhaps it needs an extra option

'LDFLAGS=-dynamic'

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Dear Satish,

I get the following when doing the ls:

[lanthale@daint103 mkl]$ ls
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/
libmkl_ao_worker.so libmkl_blacs_sgimpt_lp64.so
libmkl_intel_lp64.alibmkl_sequential.so
libmkl_avx2.so  libmkl_blas95_ilp64.a libmkl_intel_lp64.so
libmkl_tbb_thread.a
libmkl_avx512_mic.solibmkl_blas95_lp64.a libmkl_intel_thread.a
libmkl_tbb_thread.so
libmkl_avx512.solibmkl_cdft_core.a libmkl_intel_thread.so
libmkl_vml_avx2.so
libmkl_avx.so   libmkl_cdft_core.so libmkl_lapack95_ilp64.a
libmkl_vml_avx512_mic.so
libmkl_blacs_intelmpi_ilp64.a   libmkl_core.a libmkl_lapack95_lp64.a
libmkl_vml_avx512.so
libmkl_blacs_intelmpi_ilp64.so  libmkl_core.so libmkl_mc3.so
libmkl_vml_avx.so
libmkl_blacs_intelmpi_lp64.alibmkl_def.so libmkl_mc.so
libmkl_vml_cmpt.so
libmkl_blacs_intelmpi_lp64.so   libmkl_gf_ilp64.a libmkl_pgi_thread.a
libmkl_vml_def.so
libmkl_blacs_openmpi_ilp64.alibmkl_gf_ilp64.so libmkl_pgi_thread.so
libmkl_vml_mc2.so
libmkl_blacs_openmpi_ilp64.so   libmkl_gf_lp64.a libmkl_rt.so
libmkl_vml_mc3.so
libmkl_blacs_openmpi_lp64.a libmkl_gf_lp64.so libmkl_scalapack_ilp64.a
libmkl_vml_mc.so
libmkl_blacs_openmpi_lp64.solibmkl_gnu_thread.a libmkl_scalapack_ilp64.so
locale
libmkl_blacs_sgimpt_ilp64.a libmkl_gnu_thread.so libmkl_scalapack_lp64.a
libmkl_blacs_sgimpt_ilp64.solibmkl_intel_ilp64.a libmkl_scalapack_lp64.so
libmkl_blacs_sgimpt_lp64.a  libmkl_intel_ilp64.so libmkl_sequential.a


Cheers,
Samuel

On 04/24/2018 03:50 PM, Satish Balay wrote:

Executing: cc  -o /tmp/petsc-tyd_Hm/config.libraries/conftest -O
/tmp/petsc-tyd_Hm/config.libraries/conftest.o
-Wl,-rpath,/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64
-L/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64
-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lstdc++ -ldl
Possible ERROR while running linker: exit code 256
stderr:
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(mkl_semaphore.o):
In function `mkl_serv_load_inspector':
mkl_semaphore.c:(.text+0x123): warning: Using 'dlopen' in statically linked
applications requires at runtime the shared libraries from the glibc version
used for linking
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(dgetrs.o):
In function `mkl_lapack_dgetrs':
dgetrs_gen.f:(.text+0x224): undefined reference to `mkl_blas_dtrsm'
dgetrs_gen.f:(.text+0x2b1): undefined reference to `mkl_blas_dtrsm'
dgetrs_gen.f:(.text+0x2e7): undefined reference to `mkl_lapack_dlaswp'
<<<<<<<<

For some reason the compiler is picking up static version of MKL - instead
of shared libraries.

What do you have for:

ls /opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Dear Satish and Stefano,

Thank you for your answers. I believe I had initially tried to use the
option
--with-blaslapack-dir=[...], instead of specifying lib and include
directly.
But that gives me an error message:

***
   UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for
details):
---
You set a value for --with-blaslapack-dir=, but
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl cannot be used
***

I attach the new configure.log. Do you think there is something wrong with
the
mkl version that I'm trying to u

Re: [petsc-users] configuration option for PETSc on cluster

[Samuel Lanthaler]

Sorry, I'm confused: Can I add this somehow when doing the ./configure? 
Or do you mean to add it when compiling ex5f with the configuration I 
had initially?


Cheers,
Samuel


On 04/24/2018 04:03 PM, Satish Balay wrote:

Ok - I do see: libmkl_intel_lp64.so libmkl_sequential.so libmkl_core.so - these 
should
be the ones that get picked up.

I do not know why 'cc' is picking up static libraries instead of the
shared ones. Perhaps it needs an extra option

'LDFLAGS=-dynamic'

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Dear Satish,

I get the following when doing the ls:

[lanthale@daint103 mkl]$ ls
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/
libmkl_ao_worker.so libmkl_blacs_sgimpt_lp64.so
libmkl_intel_lp64.alibmkl_sequential.so
libmkl_avx2.so  libmkl_blas95_ilp64.a libmkl_intel_lp64.so
libmkl_tbb_thread.a
libmkl_avx512_mic.solibmkl_blas95_lp64.a libmkl_intel_thread.a
libmkl_tbb_thread.so
libmkl_avx512.solibmkl_cdft_core.a libmkl_intel_thread.so
libmkl_vml_avx2.so
libmkl_avx.so   libmkl_cdft_core.so libmkl_lapack95_ilp64.a
libmkl_vml_avx512_mic.so
libmkl_blacs_intelmpi_ilp64.a   libmkl_core.a libmkl_lapack95_lp64.a
libmkl_vml_avx512.so
libmkl_blacs_intelmpi_ilp64.so  libmkl_core.so libmkl_mc3.so
libmkl_vml_avx.so
libmkl_blacs_intelmpi_lp64.alibmkl_def.so libmkl_mc.so
libmkl_vml_cmpt.so
libmkl_blacs_intelmpi_lp64.so   libmkl_gf_ilp64.a libmkl_pgi_thread.a
libmkl_vml_def.so
libmkl_blacs_openmpi_ilp64.alibmkl_gf_ilp64.so libmkl_pgi_thread.so
libmkl_vml_mc2.so
libmkl_blacs_openmpi_ilp64.so   libmkl_gf_lp64.a libmkl_rt.so
libmkl_vml_mc3.so
libmkl_blacs_openmpi_lp64.a libmkl_gf_lp64.so libmkl_scalapack_ilp64.a
libmkl_vml_mc.so
libmkl_blacs_openmpi_lp64.solibmkl_gnu_thread.a libmkl_scalapack_ilp64.so
locale
libmkl_blacs_sgimpt_ilp64.a libmkl_gnu_thread.so libmkl_scalapack_lp64.a
libmkl_blacs_sgimpt_ilp64.solibmkl_intel_ilp64.a libmkl_scalapack_lp64.so
libmkl_blacs_sgimpt_lp64.a  libmkl_intel_ilp64.so libmkl_sequential.a


Cheers,
Samuel

On 04/24/2018 03:50 PM, Satish Balay wrote:

Executing: cc  -o /tmp/petsc-tyd_Hm/config.libraries/conftest -O
/tmp/petsc-tyd_Hm/config.libraries/conftest.o
-Wl,-rpath,/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64
-L/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64
-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lstdc++ -ldl
Possible ERROR while running linker: exit code 256
stderr:
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(mkl_semaphore.o):
In function `mkl_serv_load_inspector':
mkl_semaphore.c:(.text+0x123): warning: Using 'dlopen' in statically linked
applications requires at runtime the shared libraries from the glibc version
used for linking
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(dgetrs.o):
In function `mkl_lapack_dgetrs':
dgetrs_gen.f:(.text+0x224): undefined reference to `mkl_blas_dtrsm'
dgetrs_gen.f:(.text+0x2b1): undefined reference to `mkl_blas_dtrsm'
dgetrs_gen.f:(.text+0x2e7): undefined reference to `mkl_lapack_dlaswp'
<<<<<<<<

For some reason the compiler is picking up static version of MKL - instead
of shared libraries.

What do you have for:

ls /opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Dear Satish and Stefano,

Thank you for your answers. I believe I had initially tried to use the
option
--with-blaslapack-dir=[...], instead of specifying lib and include
directly.
But that gives me an error message:

***
   UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for
details):
---
You set a value for --with-blaslapack-dir=, but
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl cannot be used
***

I attach the new configure.log. Do you think there is something wrong with
the
mkl version that I'm trying to use? Or is it only related to petsc?

Cheers,
Samuel


On 04/24/2018 02:55 PM, Satish Balay wrote:

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Hi there,

I was wondering if you could help me with the correct configuration of
PETSc-dev version on a cluster (https://www.cscs.ch/)? I'm not sure which
information would be useful to you, but basically the problem seems to be
in
correctly compiling it with intel compiler and the existing mkl library.

The pre-installed mkl version is

intel/17.0.4.196

I have tried various things and finally, I got it at least to compile
with
the
following options (options chosen in reference to the mkl link
advisor...):

./configure --with-fc=ftn --with-cc=cc --with-cxx=CC --with-debugging=0
--with-scalar

Re: [petsc-users] configuration option for PETSc on cluster

[Samuel Lanthaler]

Dear Satish,

I get the following when doing the ls:

[lanthale@daint103 mkl]$ ls 
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/
libmkl_ao_worker.so libmkl_blacs_sgimpt_lp64.so 
libmkl_intel_lp64.alibmkl_sequential.so
libmkl_avx2.so  libmkl_blas95_ilp64.a 
libmkl_intel_lp64.so   libmkl_tbb_thread.a
libmkl_avx512_mic.solibmkl_blas95_lp64.a 
libmkl_intel_thread.a  libmkl_tbb_thread.so
libmkl_avx512.solibmkl_cdft_core.a 
libmkl_intel_thread.so libmkl_vml_avx2.so
libmkl_avx.so   libmkl_cdft_core.so 
libmkl_lapack95_ilp64.alibmkl_vml_avx512_mic.so
libmkl_blacs_intelmpi_ilp64.a   libmkl_core.a libmkl_lapack95_lp64.a 
libmkl_vml_avx512.so
libmkl_blacs_intelmpi_ilp64.so  libmkl_core.so 
libmkl_mc3.so  libmkl_vml_avx.so
libmkl_blacs_intelmpi_lp64.alibmkl_def.so libmkl_mc.so   
libmkl_vml_cmpt.so
libmkl_blacs_intelmpi_lp64.so   libmkl_gf_ilp64.a 
libmkl_pgi_thread.alibmkl_vml_def.so
libmkl_blacs_openmpi_ilp64.alibmkl_gf_ilp64.so 
libmkl_pgi_thread.so   libmkl_vml_mc2.so
libmkl_blacs_openmpi_ilp64.so   libmkl_gf_lp64.a 
libmkl_rt.so   libmkl_vml_mc3.so
libmkl_blacs_openmpi_lp64.a libmkl_gf_lp64.so 
libmkl_scalapack_ilp64.a   libmkl_vml_mc.so
libmkl_blacs_openmpi_lp64.solibmkl_gnu_thread.a 
libmkl_scalapack_ilp64.so  locale

libmkl_blacs_sgimpt_ilp64.a libmkl_gnu_thread.so libmkl_scalapack_lp64.a
libmkl_blacs_sgimpt_ilp64.solibmkl_intel_ilp64.a 
libmkl_scalapack_lp64.so

libmkl_blacs_sgimpt_lp64.a  libmkl_intel_ilp64.so libmkl_sequential.a


Cheers,
Samuel

On 04/24/2018 03:50 PM, Satish Balay wrote:

Executing: cc  -o /tmp/petsc-tyd_Hm/config.libraries/conftest -O 
/tmp/petsc-tyd_Hm/config.libraries/conftest.o  
-Wl,-rpath,/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64 
-L/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64 
-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lstdc++ -ldl
Possible ERROR while running linker: exit code 256
stderr:
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(mkl_semaphore.o):
 In function `mkl_serv_load_inspector':
mkl_semaphore.c:(.text+0x123): warning: Using 'dlopen' in statically linked 
applications requires at runtime the shared libraries from the glibc version 
used for linking
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/libmkl_core.a(dgetrs.o):
 In function `mkl_lapack_dgetrs':
dgetrs_gen.f:(.text+0x224): undefined reference to `mkl_blas_dtrsm'
dgetrs_gen.f:(.text+0x2b1): undefined reference to `mkl_blas_dtrsm'
dgetrs_gen.f:(.text+0x2e7): undefined reference to `mkl_lapack_dlaswp'
<<<<<<<<

For some reason the compiler is picking up static version of MKL - instead of 
shared libraries.

What do you have for:

ls /opt/intel/compilers_and_libraries_2017.4.196/linux/mkl/lib/intel64/

Satish

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Dear Satish and Stefano,

Thank you for your answers. I believe I had initially tried to use the option
--with-blaslapack-dir=[...], instead of specifying lib and include directly.
But that gives me an error message:

***
  UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for
details):
---
You set a value for --with-blaslapack-dir=, but
/opt/intel/compilers_and_libraries_2017.4.196/linux/mkl cannot be used
***

I attach the new configure.log. Do you think there is something wrong with the
mkl version that I'm trying to use? Or is it only related to petsc?

Cheers,
Samuel


On 04/24/2018 02:55 PM, Satish Balay wrote:

On Tue, 24 Apr 2018, Samuel Lanthaler wrote:


Hi there,

I was wondering if you could help me with the correct configuration of
PETSc-dev version on a cluster (https://www.cscs.ch/)? I'm not sure which
information would be useful to you, but basically the problem seems to be
in
correctly compiling it with intel compiler and the existing mkl library.

The pre-installed mkl version is

intel/17.0.4.196

I have tried various things and finally, I got it at least to compile with
the
following options (options chosen in reference to the mkl link advisor...):

./configure --with-fc=ftn --with-cc=cc --with-cxx=CC --with-debugging=0
--with-scalar-type=complex --download-scalapack --download-mumps=yes
--download-superlu --with-batch --known-mpi-shared-libraries=1
*--with-blaslapack-lib=*" ${MKLROOT}/lib/intel64/libmkl_blas95_lp64.a
${MKLROOT}/lib/intel64/libmkl_lapack95_lp64.a -Wl,--start-group
${MKLROOT}/lib/intel64/libmkl_intel_lp64.a
${MKLROOT}/lib/intel64/libmkl_sequential.a
${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl"
*--wit

Re: [petsc-users] VecGetArrayReadF90 with gfortran

[Samuel Lanthaler]

Dear Satish,

Yes, I used the makefile you provided.


On 23.03.2018 19:00, Satish Balay wrote:

Have you compiled the my modified code with petsc makefile that I provided?

I've check with - but I can't reproduce this problem.

balay@es^/sandbox/balay/test $ gfortran --version
GNU Fortran (GCC) 6.2.0
Copyright (C) 2016 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Satish


On Fri, 23 Mar 2018, Samuel Lanthaler wrote:


Thank you very much for your effort Satish. It appears that even in your
version of the example, I'm having the same problem with a pointer that is not
properly assigned. Seems weird, I also tried recompiling petsc, but that gave
me the same. The configure file for petsc is also attached.

Btw. It appears the latest version of gfortran I have access to on my work
station is gcc-6.3.


mpif90 --version

GNU Fortran (GCC) 6.3.0
Copyright (C) 2016 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Does anyone maybe have any ideas?

Thanks,

Samuel


On 03/23/2018 05:42 PM, Satish Balay wrote:

Attaching the modified version that works with complex [using complex build
of petsc]

BTW: I'm using the following gfortran [but I don't think the version makes
any difference] with petsc-3.8 (maint branch):

$ gfortran --version
GNU Fortran (GCC) 8.0.1 20180222 (Red Hat 8.0.1-0.16)

Satish



balay@asterix /home/balay/download-pine
$ make ex2f
mpif90 -c -Wall -ffree-line-length-0 -Wno-unused-dummy-argument -g
-I/home/balay/petsc/include -I/home/balay/petsc/arch-cmplx/include-o
ex2f.o ex2f.F90
mpif90 -Wall -ffree-line-length-0 -Wno-unused-dummy-argument -g   -o ex2f
ex2f.o  -Wl,-rpath,/home/balay/petsc/arch-cmplx/lib
-L/home/balay/petsc/arch-cmplx/lib
-Wl,-rpath,/home/balay/soft/mpich-3.3b1/lib
-L/home/balay/soft/mpich-3.3b1/lib
-Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/8
-L/usr/lib/gcc/x86_64-redhat-linux/8 -lpetsc -llapack -lblas -lX11 -lpthread
-lm -lmpifort -lgfortran -lm -lgfortran -lm -lquadmath -lmpicxx -lstdc++ -lm
-Wl,-rpath,/home/balay/soft/mpich-3.3b1/lib
-L/home/balay/soft/mpich-3.3b1/lib
-Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/8
-L/usr/lib/gcc/x86_64-redhat-linux/8 -ldl
-Wl,-rpath,/home/balay/soft/mpich-3.3b1/lib -lmpi -lgcc_s -ldl
/usr/bin/rm -f ex2f.o
balay@asterix /home/balay/download-pine
$ ./ex2f
   setting values
   assembling vector
   copying vector to array
   N =   10
   m =   10
   getting vecX size
   getting ownership range
   ownership range:0  10
   copy values to global array
   vecX_pt:
   SIZE =   10
(1.,0.)
(2.,0.)
(3.,0.)
(4.,0.)
(5.,0.)
(6.,0.)
(7.,0.)
(8.,0.)
(9.,0.)
(10.000,0.)
   copy_buffer =
(1.,0.)
(2.,0.)
(3.,0.)
(4.,0.)
(5.,0.)
(6.,0.)
(7.,0.)
(8.,0.)
(9.,0.)
(10.000,0.)
   calling allreduce
   buffer =
(1.,0.)
(2.,0.)
(3.,0.)
(4.,0.)
(5.,0.)
(6.,0.)
(7.,0.)
(8.,0.)
(9.,0.)
(10.000,0.)
balay@asterix /home/balay/download-pine
$

On Fri, 23 Mar 2018, Satish Balay wrote:


I'm attaching the modified code - that appears to work fine for me. My
changes are:

1. use .F90 suffix - so the compiler does the f90 freeform/preprocessing
automatically
2. remove slepc references and use petsc only
3. compile with a petsc formatted

Re: [petsc-users] VecGetArrayReadF90 with gfortran

[Samuel Lanthaler]

Thanks for your answers.

@Timothee: I think normally, it should be fine to use the pointer like 
an array. At least, I have done the same in some other cases. But I can 
try anyways.


@Barry: Replacing the copy by a for-loop does not seem to help in this 
case: the output pointer has LBOUND=1 and UBOUND=0, SIZE=0, so that 
trying to access any index leads to a segfault.



On 03/23/2018 04:30 PM, Smith, Barry F. wrote:

Use a loop to do the copy

At line 127 of file allreduce.f90
Fortran runtime error: Array bound mismatch for dimension 1 of array 
'copy_buffer' (10/0)




On Mar 23, 2018, at 9:10 AM, Samuel Lanthaler <s.lantha...@gmail.com> wrote:

Hi all,

I am having trouble using the function VecGetArrayReadF90 under gfortran. I 
have created a minimal example and put it in the attachment. Basically, it 
appears that the input pointer is empty upon return from VecGetArrayReadF90. 
The code runs fine, when compiled with ifort, so I don't know what's going on.

I am attaching the code and makefile that I used, as well as the command line 
output that I can see. Though the errors are just to do with the fact that the 
output pointer is empty, so the command line output is not all that helpful...

Thanks in advance for your help!
Sam

[petsc-users] VecGetArrayReadF90 with gfortran

[Samuel Lanthaler]

Hi all,

I am having trouble using the function VecGetArrayReadF90 under gfortran. I
have created a minimal example and put it in the attachment. Basically, it
appears that the input pointer is empty upon return from
VecGetArrayReadF90. The code runs fine, when compiled with ifort, so I
don't know what's going on.

I am attaching the code and makefile that I used, as well as the command
line output that I can see. Though the errors are just to do with the fact
that the output pointer is empty, so the command line output is not all
that helpful...

Thanks in advance for your help!
Sam
 setting values
 assembling vector
 copying vector to array
 N =   10
 m =   10
 getting vecX size
 getting ownership range
 ownership range:0  10
 copy values to global array
 vecX_pt:
 SIZE =0
At line 127 of file allreduce.f90
Fortran runtime error: Array bound mismatch for dimension 1 of array 
'copy_buffer' (10/0)

Error termination. Backtrace:
#0  0x402621 in petscvec2array_
at /home/lanthale/Progs/test_PetscVec2Array/allreduce.f90:127
#1  0x403550 in allreduce
at /home/lanthale/Progs/test_PetscVec2Array/allreduce.f90:58
#2  0x4035d8 in main
at /home/lanthale/Progs/test_PetscVec2Array/allreduce.f90:4


makefile
Description: Binary data
PROGRAM allreduce
#include "slepc/finclude/slepc.h"
  ! --
  USE slepcsys
  USE slepceps
  IMPLICIT NONE
  
  ! --- pure PETSc
  PetscErrorCode :: ierr
  PetscScalar, ALLOCATABLE :: vals(:)
  PetscInt, ALLOCATABLE :: idxn(:)
  PetscInt :: m,i
  ! --- SLEPc
  Vec :: vecX! eigenvector (real/imaginary parts)
  !
  COMPLEX, ALLOCATABLE :: vecXcopy(:)
  INTEGER :: N

  ! initialize SLEPc
  CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)

  ! Set up a new matrix
  m = 10
  ALLOCATE(vals(10),idxn(10))

  !
  CALL VecCreate(PETSC_COMM_WORLD,vecX,ierr); 
  CALL VecSetType(vecX,VECMPI,ierr); 
  CALL VecSetSizes(vecX,PETSC_DECIDE,m,ierr); 

  ! set values of vector
  DO i=1,m
 vals(i) = i
 idxn(i) = i-1
  END DO
  !
  PRINT*,'setting values'
  CALL VecSetValues(vecX,m,idxn,vals,INSERT_VALUES,ierr); 

  ! assemble vector
  PRINT*,'assembling vector'
  CALL VecAssemblyBegin(vecX,ierr); 
  CALL VecAssemblyEnd(vecX,ierr); 

  ! 
  PRINT*,'copying vector to array'
  N = m
  PRINT*,'N = ',N
  PRINT*,'m = ',m
  ALLOCATE(vecXcopy(N))
  CALL PetscVec2Array(vecX,vecXcopy,N)

  !
  CALL VecDestroy(vecX,ierr); 
  CALL SlepcFinalize(ierr); 


END PROGRAM allreduce


  !> Copy Petsc vector to fortran array
  SUBROUTINE PetscVec2Array(vecX,buffer,N)
#include "slepc/finclude/slepc.h"
! --
!USE petscsys
!USE petscvec
! --
USE slepcsys
USE slepceps
IMPLICIT NONE
Vec, INTENT(in) :: vecX
COMPLEX, INTENT(inout) :: buffer(0:N-1)
INTEGER, INTENT(in) :: N
!
INTEGER :: istart, iend, Np, i,j
INTEGER :: me2,nprocs
COMPLEX, ALLOCATABLE :: copy_buffer(:), copy_buffer2(:)
!
PetscScalar, POINTER :: vecX_pt(:) => NULL()
PetscInt :: petsc_istart, petsc_iend, petsc_N
PetscErrorCode :: ierr

!
CALL MPI_COMM_RANK(MPI_COMM_WORLD,me2,ierr);
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,Nprocs,ierr);

! initailize buffer to 0.
buffer = 0.
ALLOCATE(copy_buffer(0:N-1))
   

! check length of petsc vector is N
PRINT*,'getting vecX size'
CALL VecGetSize(vecX,petsc_N,ierr); CHKERRQ(ierr)
Np = petsc_N
IF(N.NE.Np) THEN
   PRINT*,'ERROR: In PetscVec2Array:'
   PRINT*,'petsc_N = ',petsc_N, Np
   PRINT*,'N   = ',N
   CALL ABORT
END IF

! vecX_pt will contain the element in the local range [istart,iend-1]
PRINT*,'getting ownership range'
CALL VecGetOwnershipRange(vecX,istart,iend,ierr); CHKERRQ(ierr)
PRINT*,'ownership range: ',istart,iend

! initialize the pointer
CALL VecGetArrayReadF90(vecX,vecX_pt,ierr); CHKERRQ(ierr)
! copy values to global array
PRINT*,'copy values to global array'

DO i=1,Nprocs
   !
   CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
   !
   IF(me2.EQ.i-1) THEN
  PRINT*,'vecX_pt:'
  PRINT*,'SIZE = ',SIZE(VecX_pt)
  DO j=1,SIZE(vecX_pt)
 PRINT*,'  ',vecX_pt(j)
  END DO
   END IF
END DO


copy_buffer(istart:iend-1) = vecX_pt(:)
! free pointer
CALL VecRestoreArrayReadF90(vecX,vecX_pt,ierr); CHKERRQ(ierr)

CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

DO i=1,Nprocs
   !
   CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
   !
   IF(me2.EQ.i-1) THEN
  PRINT*,'copy_buffer = '
  DO j=istart,iend-1
 PRINT*,'  ',copy_buffer(j)
  END DO
   END IF
END DO

CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

! combine all components from all processes
PRINT*,'calling allreduce'
CALL MPI_ALLREDUCE(copy_buffer,buffer,N,MPI_DOUBLE_COMPLEX,MPI_SUM,MPI_COMM_WORLD,ierr)
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

DO i=1,Nprocs
   !
   CALL 

Re: [petsc-users] make: *** No rule to make target `chk_makej', needed by `all'. Stop.

[Samuel Lanthaler]

Ah, now I undestood what you meant... sorry. Thank you, Jose!


On 08.03.2018 22:26, Jed Brown wrote:

Thanks, Jose.

"Jose E. Roman" <jro...@dsic.upv.es> writes:


Update the SLEPc repository and try again.
When a change in petsc-master affects SLEPc, you have to allow one or two days 
for us to re-synchronize.
Jose


El 9 mar 2018, a las 3:59, Samuel Lanthaler <s.lantha...@gmail.com> escribió:

Hi,

I've just been trying to compile SLEPC on a cluster, and am getting the following error 
message after I type "make":


./configure

Checking environment... done
Checking PETSc installation... done
Checking LAPACK library... done
Writing various configuration files... done
Generating Fortran stubs... done

===
SLEPc Configuration
===

SLEPc directory:
/marconi/home/userexternal/slanthal/slepc-git
  It is a git repository on branch: master
PETSc directory:
/marconi/home/userexternal/slanthal/petsc-git
  It is a git repository on branch: master
Architecture "arch-complex" with double precision complex numbers

xxx=xxx
Configure stage complete. Now build the SLEPc library with (gnumake build):
   make SLEPC_DIR=$PWD PETSC_DIR=/marconi/home/userexternal/slanthal/petsc-git 
PETSC_ARCH=arch-complex
xxx=xxx


make SLEPC_DIR=$PWD PETSC_DIR=/marconi/home/userexternal/slanthal/petsc-git 
PETSC_ARCH=arch-complex

make: *** No rule to make target `chk_makej', needed by `all'. Stop.

Does someone maybe know what the cause of this could be? Make is correct in 
that there is no rule for chk_makej, but I'm assuming it should really be some 
kind of function? However, since I don't understand what the dependency on 
chk_makej exactly does, I'm a bit lost...

Thanks for you help in advance!
Samuel

Re: [petsc-users] make: *** No rule to make target `chk_makej', needed by `all'. Stop.

[Samuel Lanthaler]

Thanks, Jose. I actually tried that before, but it turns out there 
probably was some kind of an issue already with the compilation of 
petsc; even though it looked like it had worked. Trying again, now.


Thanks,

Samuel


On 08.03.2018 22:14, Jose E. Roman wrote:

Update the SLEPc repository and try again.
When a change in petsc-master affects SLEPc, you have to allow one or two days 
for us to re-synchronize.
Jose


El 9 mar 2018, a las 3:59, Samuel Lanthaler <s.lantha...@gmail.com> escribió:

Hi,

I've just been trying to compile SLEPC on a cluster, and am getting the following error 
message after I type "make":


./configure

Checking environment... done
Checking PETSc installation... done
Checking LAPACK library... done
Writing various configuration files... done
Generating Fortran stubs... done

===
SLEPc Configuration
===

SLEPc directory:
/marconi/home/userexternal/slanthal/slepc-git
  It is a git repository on branch: master
PETSc directory:
/marconi/home/userexternal/slanthal/petsc-git
  It is a git repository on branch: master
Architecture "arch-complex" with double precision complex numbers

xxx=xxx
Configure stage complete. Now build the SLEPc library with (gnumake build):
   make SLEPC_DIR=$PWD PETSC_DIR=/marconi/home/userexternal/slanthal/petsc-git 
PETSC_ARCH=arch-complex
xxx=xxx


make SLEPC_DIR=$PWD PETSC_DIR=/marconi/home/userexternal/slanthal/petsc-git 
PETSC_ARCH=arch-complex

make: *** No rule to make target `chk_makej', needed by `all'. Stop.

Does someone maybe know what the cause of this could be? Make is correct in 
that there is no rule for chk_makej, but I'm assuming it should really be some 
kind of function? However, since I don't understand what the dependency on 
chk_makej exactly does, I'm a bit lost...

Thanks for you help in advance!
Samuel

[petsc-users] make: *** No rule to make target `chk_makej', needed by `all'. Stop.

[Samuel Lanthaler]

Hi,

I've just been trying to compile SLEPC on a cluster, and am getting the 
following error message after I type "make":


> ./configure
Checking environment... done
Checking PETSc installation... done
Checking LAPACK library... done
Writing various configuration files... done
Generating Fortran stubs... done

===
SLEPc Configuration
===

SLEPc directory:
 /marconi/home/userexternal/slanthal/slepc-git
  It is a git repository on branch: master
PETSc directory:
 /marconi/home/userexternal/slanthal/petsc-git
  It is a git repository on branch: master
Architecture "arch-complex" with double precision complex numbers

xxx=xxx
 Configure stage complete. Now build the SLEPc library with (gnumake 
build):
   make SLEPC_DIR=$PWD 
PETSC_DIR=/marconi/home/userexternal/slanthal/petsc-git 
PETSC_ARCH=arch-complex

xxx=xxx

> make SLEPC_DIR=$PWD 
PETSC_DIR=/marconi/home/userexternal/slanthal/petsc-git 
PETSC_ARCH=arch-complex

make: *** No rule to make target `chk_makej', needed by `all'. Stop.

Does someone maybe know what the cause of this could be? Make is correct 
in that there is no rule for chk_makej, but I'm assuming it should 
really be some kind of function? However, since I don't understand what 
the dependency on chk_makej exactly does, I'm a bit lost...


Thanks for you help in advance!
Samuel

Re: [petsc-users] Using SLEPC: Fortran call to NEPNLEIGSSetSingularitiesFunction (similar to SNESSetFunction?)

[Samuel Lanthaler]

Amazing, thank you!
Samuel


On 01/11/2018 02:48 PM, Jose E. Roman wrote:

The Fortran stub is now in master. I have also completed the example.
https://bitbucket.org/slepc/slepc/commits/acaefe873
Jose


El 10 ene 2018, a las 17:29, Samuel Lanthaler <s.lantha...@gmail.com> escribió:

Ah, I see.
Great, thank you, Jose!

Cheers,
Samuel


On 01/10/2018 05:23 PM, Jose E. Roman wrote:

This function is not implemented in Fortran. I will add it and let you know.
Jose



El 10 ene 2018, a las 17:11, Samuel Lanthaler <s.lantha...@gmail.com> escribió:

Hi there,

I'm stuck getting a call to the SLEPC routine 
"NEPNLEIGSSetSingularitiesFunction" to work from within Fortran.
To be more precise: To get started with the solution of a non-linear 
eigenvalue-problem, I am first trying to write a Fortran version of 
slepc-3.8.1/.../nep/examples/tutorials/ex27.c, and then use it as a template 
for my own code. But I haven't gotten far...

Let me attach the source of the current state of my attempt to translate ex27.c 
to Fortran code. When trying to compile the attached source, my compiler 
complains:

nleigs.o: In function `MAIN__':
nleigs.f90:(.text+0x213): undefined reference to 
`nepnleigssetsingularitiesfunction_'

Does someone maybe see what I'm doing wrong? Do I maybe need to add a USE 
statement in addition to

USE slepcsys
USE slepcnep

?
I'm quite confused about the C++ pointers and how to deal with them when 
calling C++ functions from Fortran, so I don't really understand what's going 
on here. I will greatly appreciate your help!

Thanks,
Samuel

Re: [petsc-users] Using SLEPC: Fortran call to NEPNLEIGSSetSingularitiesFunction (similar to SNESSetFunction?)

[Samuel Lanthaler]

Ah, I see.
Great, thank you, Jose!

Cheers,
Samuel


On 01/10/2018 05:23 PM, Jose E. Roman wrote:

This function is not implemented in Fortran. I will add it and let you know.
Jose



El 10 ene 2018, a las 17:11, Samuel Lanthaler <s.lantha...@gmail.com> escribió:

Hi there,

I'm stuck getting a call to the SLEPC routine 
"NEPNLEIGSSetSingularitiesFunction" to work from within Fortran.
To be more precise: To get started with the solution of a non-linear 
eigenvalue-problem, I am first trying to write a Fortran version of 
slepc-3.8.1/.../nep/examples/tutorials/ex27.c, and then use it as a template 
for my own code. But I haven't gotten far...

Let me attach the source of the current state of my attempt to translate ex27.c 
to Fortran code. When trying to compile the attached source, my compiler 
complains:

nleigs.o: In function `MAIN__':
nleigs.f90:(.text+0x213): undefined reference to 
`nepnleigssetsingularitiesfunction_'

Does someone maybe see what I'm doing wrong? Do I maybe need to add a USE 
statement in addition to

USE slepcsys
USE slepcnep

?
I'm quite confused about the C++ pointers and how to deal with them when 
calling C++ functions from Fortran, so I don't really understand what's going 
on here. I will greatly appreciate your help!

Thanks,
Samuel

[petsc-users] Using SLEPC: Fortran call to NEPNLEIGSSetSingularitiesFunction (similar to SNESSetFunction?)

[Samuel Lanthaler]

Hi there,

I'm stuck getting a call to the SLEPC routine 
"NEPNLEIGSSetSingularitiesFunction" to work from within Fortran.
To be more precise: To get started with the solution of a non-linear 
eigenvalue-problem, I am first trying to write a Fortran version of 
slepc-3.8.1/.../nep/examples/tutorials/ex27.c, and then use it as a 
template for my own code. But I haven't gotten far...


Let me attach the source of the current state of my attempt to translate 
ex27.c to Fortran code. When trying to compile the attached source, my 
compiler complains:


nleigs.o: In function `MAIN__':
nleigs.f90:(.text+0x213): undefined reference to 
`nepnleigssetsingularitiesfunction_'


Does someone maybe see what I'm doing wrong? Do I maybe need to add a 
USE statement in addition to


USE slepcsys
USE slepcnep

?
I'm quite confused about the C++ pointers and how to deal with them when 
calling C++ functions from Fortran, so I don't really understand what's 
going on here. I will greatly appreciate your help!


Thanks,
Samuel
! 
!user-defined context
! 
MODULE solver_context
#include "slepc/finclude/slepc.h"
  ! --
  USE slepcsys
  USE slepcnep
  ! --
  IMPLICIT NONE

  TYPE :: MatCtx
 PetscScalar :: lambda,kappa
 PetscReal :: h
  END TYPE MatCtx

END MODULE solver_context

! 
!main program
! 
!
!   Solve T(lambda)x=0 using NLEIGS solver
!  with T(lambda) = -D+sqrt(lambda)*I
!  where D is the Laplacian operator in 1 dimension
!  and with the interpolation interval [.01,16]
!
PROGRAM main
#include "slepc/finclude/slepc.h"
  ! --
  USE slepcsys
  USE slepcnep
  ! --
  USE solver_context
  IMPLICIT NONE
  NEP :: nep
  Mat :: F,A(2)
  NEPType :: type
  PetscInt :: n=100,nev,Istart,Iend,i
  PetscErrorCode :: ierr
  PetscBool :: terse,split=PETSC_TRUE
  RG :: rg
  FN :: fn(2)
  PetscScalar :: coeffs
  PetscBool :: flg
  ! 
  CHARACTER(LEN=128) :: string
  TYPE(MatCtx),TARGET :: ctx
  TYPE(MatCtx),POINTER :: ctx_pt

  ! NOTE: Any user-defined Fortran routines (such as ComputeSingularities)
  !   MUST be declared as external.
  external ComputeSingularities

  !
  CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
  CALL PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-n",n,flg,ierr)
  CALL PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-split",split,flg,ierr)
  WRITE(string,*) 'Square root eigenproblem, n=',n,'\n In split form: ',split,'\n'
  CALL PetscPrintf(PETSC_COMM_WORLD,TRIM(string),ierr)
  
  ! --
  ! Create nonlinear eigensolver context
  ! --
  CALL NEPCreate(PETSC_COMM_WORLD,nep,ierr)

  ! --
  ! Select the NLEIGS solver and set required options for it
  ! --
  CALL NEPSetType(nep,NEPNLEIGS,ierr)
  ctx_pt => ctx
  CALL NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,ctx_pt,ierr)


  !
  CALL SlepcFinalize(ierr)

END PROGRAM main


! ComputeSingularities --
! Computes maxnp points (at most) in the complex plane where the
! function T(.) is not analytic.
SUBROUTINE ComputeSingularities(nep,maxnp,xi,ctx_pt,ierr)
#include "slepc/finclude/slepc.h"
  USE solver_context
  ! ---
  USE slepcsys
  USE slepcnep
  IMPLICIT NONE
  NEP :: nep
  PetscInt :: maxnp
  PetscScalar :: xi(0:maxnp-1)
  TYPE(MatCtx), POINTER :: ctx_pt
  PetscErrorCode :: ierr
  ! ---
  PetscReal :: h
  PetscInt :: i

  h = 11.0/(maxnp-1)
  xi(0) = -1e-5
  xi(maxnp-1) = -1e+6
  DO i=1,maxnp-2
 xi(i) = -10**(-5+h*i)
  END DO

END SUBROUTINE ComputeSingularities

Re: [petsc-users] MatCreateShell, MatShellGetContext, MatShellSetContext in fortran

[Samuel Lanthaler]

Ah, silly me... Of course, if the program can't actually see the 
interface then it makes sense that it won't work.

As always, thanks a lot for your help, Barry!

Samuel


On 12/13/2017 04:01 AM, Smith, Barry F. wrote:


  Samuel,

I have attached a fixed up version of your example that works 
correctly. Your basic problem was that you did not "use" the final 
module in the main program so it did know the interfaces for the 
routines.   I have included this example in PETSc for others to 
benefit from in the branch barry/add-f90-matshellgetcontext-example


   Thanks for the question,

   Barry




> On Dec 12, 2017, at 3:52 AM, Samuel Lanthaler 
<s.lantha...@gmail.com> wrote:

>
> Let me also add a minimal example (relying only on petsc), which 
leads to the same error message on my machine. Again, what I'm trying 
to do is very simple:

>• MatCreateShell => Initialize Mat :: F
>• MatShellSetContext => set the context of F to ctxF (looking 
at the petsc source code, this call actually seems to be superfluous, 
but nevermind)
>• MatShellGetContext => get the pointer ctxF_pt to point to 
the matrix context

> I'm getting an error message in the third step.
>
> [0]PETSC ERROR: - Error Message 
--

> [0]PETSC ERROR: Null argument, when expecting valid pointer
> [0]PETSC ERROR: Null Pointer: Parameter # 2
> Again, thanks for your help!
> Cheers,
> Samuel
>
> On 12/11/2017 07:41 PM, Samuel Lanthaler wrote:
>> Dear petsc-/slepc-users,
>>
>> I have been trying to understand matrix-free/shell matrices in 
PETSc for eventual use in solving a non-linear eigenvalue problem 
using SLEPC. But I seem to be having trouble with calls to 
MatShellGetContext. As far as I understand, this function should 
initialize a pointer (second argument) so that the subroutine output 
will point to the context associated with my shell-matrix (let's say 
of TYPE(MatCtx))? When calling that subroutine, I get the following 
error message:

>>
>> [0]PETSC ERROR: - Error Message 
--

>> [0]PETSC ERROR: Null argument, when expecting valid pointer
>> [0]PETSC ERROR: Null Pointer: Parameter # 2
>>
>> In my code, the second input argument to the routine is a 
null-pointer of TYPE(MatCtx),POINTER :: arg2. Which the error message 
appears to be unhappy with. I noticed that there is no error message 
if I instead pass an object TYPE(MatCtx) :: arg2 to the routine... 
which doesn't really make sense to me? Could someone maybe explain to 
me what is going on, here?

>>
>> Just in case, let me also attach my concrete example code (it is 
supposed to be a Fortran version of the slepc-example in 
slepc-3.8.1/src/nep/examples/tutorials/ex21.c). I have added an extra 
call to MatShellGetContext on line 138, after the function and 
jacobian should supposedly have been set up.

>>
>> Thanks a lot for your help!
>>
>> Cheers,
>>
>> Samuel
>>
>
> 

Re: [petsc-users] MatCreateShell, MatShellGetContext, MatShellSetContext in fortran

[Samuel Lanthaler]

Let me also add a minimal example (relying only on petsc), which leads 
to the same error message on my machine. Again, what I'm trying to do is 
very simple:


1. MatCreateShell => Initialize Mat :: F
2. MatShellSetContext => set the context of F to ctxF (looking at the
   petsc source code, this call actually seems to be superfluous, but
   nevermind)
3. MatShellGetContext => get the pointer ctxF_pt to point to the matrix
   context

I'm getting an error message in the third step.

[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Pointer: Parameter # 2

Again, thanks for your help!

Cheers,
Samuel


On 12/11/2017 07:41 PM, Samuel Lanthaler wrote:

Dear petsc-/slepc-users,

I have been trying to understand matrix-free/shell matrices in PETSc 
for eventual use in solving a non-linear eigenvalue problem using 
SLEPC. But I seem to be having trouble with calls to 
MatShellGetContext. As far as I understand, this function should 
initialize a pointer (second argument) so that the subroutine output 
will point to the context associated with my shell-matrix (let's say 
of TYPE(MatCtx))? When calling that subroutine, I get the following 
error message:


[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Pointer: Parameter # 2

In my code, the second input argument to the routine is a null-pointer 
of TYPE(MatCtx),POINTER :: arg2. Which the error message appears to be 
unhappy with. I noticed that there is no error message if I instead 
pass an object TYPE(MatCtx) :: arg2 to the routine... which doesn't 
really make sense to me? Could someone maybe explain to me what is 
going on, here?


Just in case, let me also attach my concrete example code (it is 
supposed to be a Fortran version of the slepc-example in 
slepc-3.8.1/src/nep/examples/tutorials/ex21.c). I have added an extra 
call to MatShellGetContext on line 138, after the function and 
jacobian should supposedly have been set up.


Thanks a lot for your help!

Cheers,

Samuel



! 
!user-defined context
! 
MODULE solver_context
#include "petsc/finclude/petsc.h"
  ! --
  USE petscsys
  USE petscmat
  ! --
  IMPLICIT NONE

  TYPE :: MatCtx
 PetscScalar :: lambda,kappa
 PetscReal :: h
  END TYPE MatCtx

END MODULE solver_context

MODULE solver_context_interfaces
  USE solver_context
  IMPLICIT NONE

  ! 
  INTERFACE MatCreateShell
 ! --
 SUBROUTINE MatCreateShell(comm,mloc,nloc,m,n,ctx,mat,ierr)
   USE solver_context
   TYPE(MPI_COMM) :: comm
   PetscInt :: mloc,nloc,m,n
   TYPE(MatCtx) :: ctx
   Mat :: mat
   PetscErrorCode :: ierr
 END SUBROUTINE MatCreateShell
  END INTERFACE MatCreateShell
  ! 

  ! 
  INTERFACE MatShellSetContext
 ! --
 SUBROUTINE MatShellSetContext(mat,ctx,ierr)
   USE solver_context
   Mat :: mat
   TYPE(MatCtx) :: ctx
   PetscErrorCode :: ierr
 END SUBROUTINE MatShellSetContext
  END INTERFACE MatShellSetContext
  ! 

  ! 
  INTERFACE MatShellGetContext
 ! --
 SUBROUTINE MatShellGetContext(mat,ctx,ierr)
   USE solver_context
   Mat :: mat
   TYPE(MatCtx), POINTER :: ctx
   PetscErrorCode :: ierr
 END SUBROUTINE MatShellGetContext
  END INTERFACE MatShellGetContext
  ! 

END MODULE solver_context_interfaces

! 
!main program
! 
PROGRAM main
#include "petsc/finclude/petsc.h"
  ! --
  USE petscsys
  USE petscmat
  ! --
  USE solver_context
  IMPLICIT NONE
  Mat :: F
  TYPE(MatCtx) :: ctxF
  TYPE(MatCtx),POINTER :: ctxF_pt
  PetscErrorCode :: ierr
  PetscInt :: n=128
  
  !
  CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)
  
  !
  ctxF%lambda = 0.0d0
  CALL MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,ctxF,F,ierr)
  CALL MatShellSetContext(F,ctxF,ierr)
  PRINT*,'ctxF%lambda = ',ctxF%lambda

  CALL MatShellGetContext(F,ctxF_pt,ierr)
  PRINT*,'ctxF_pt%lambda = ',ctxF_pt%lambda 

  !
  CALL PetscFinalize(ierr)

END PROGRAM main

[petsc-users] MatCreateShell, MatShellGetContext, MatShellSetContext in fortran

[Samuel Lanthaler]

Dear petsc-/slepc-users,

I have been trying to understand matrix-free/shell matrices in PETSc for 
eventual use in solving a non-linear eigenvalue problem using SLEPC. But 
I seem to be having trouble with calls to MatShellGetContext. As far as 
I understand, this function should initialize a pointer (second 
argument) so that the subroutine output will point to the context 
associated with my shell-matrix (let's say of TYPE(MatCtx))? When 
calling that subroutine, I get the following error message:


[0]PETSC ERROR: - Error Message 
--

[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Pointer: Parameter # 2

In my code, the second input argument to the routine is a null-pointer 
of TYPE(MatCtx),POINTER :: arg2. Which the error message appears to be 
unhappy with. I noticed that there is no error message if I instead pass 
an object TYPE(MatCtx) :: arg2 to the routine... which doesn't really 
make sense to me? Could someone maybe explain to me what is going on, here?


Just in case, let me also attach my concrete example code (it is 
supposed to be a Fortran version of the slepc-example in 
slepc-3.8.1/src/nep/examples/tutorials/ex21.c). I have added an extra 
call to MatShellGetContext on line 138, after the function and jacobian 
should supposedly have been set up.


Thanks a lot for your help!

Cheers,

Samuel

! 
!user-defined context
! 
MODULE solver_context
#include "slepc/finclude/slepc.h"
  ! --
  USE slepcsys
  USE slepceps
  USE slepcnep
  ! --
  IMPLICIT NONE
  EXTERNAL :: MatMult_Fun, MatGetDiagonal_Fun, MatDestroy_Fun, MatDuplicate_Fun
  EXTERNAL :: MatMult_Jac, MatGetDiagonal_Jac, MatDestroy_Jac

  TYPE :: ApplicationCtx
 PetscScalar :: kappa
 PetscReal :: h
  END TYPE ApplicationCtx


  TYPE :: MatCtx
 PetscScalar :: lambda,kappa
 PetscReal :: h
  END TYPE MatCtx

END MODULE solver_context

MODULE solver_context_interfaces
  USE solver_context
  IMPLICIT NONE

  ! 
  INTERFACE MatShellSetContext
 ! --
 SUBROUTINE MatShellSetContext(mat,ctx,ierr)
   USE solver_context
   Mat :: mat
   TYPE(MatCtx) :: ctx
   PetscErrorCode :: ierr
 END SUBROUTINE MatShellSetContext
  END INTERFACE MatShellSetContext
  ! 

  ! 
  INTERFACE MatShellGetContext
 ! --
 SUBROUTINE MatShellGetContext(mat,ctx,ierr)
   USE solver_context
   Mat :: mat
   TYPE(MatCtx), POINTER :: ctx
   PetscErrorCode :: ierr
 END SUBROUTINE MatShellGetContext
  END INTERFACE MatShellGetContext
  ! 

END MODULE solver_context_interfaces

! 
!main program
! 
PROGRAM main
#include "slepc/finclude/slepc.h"
  ! --
  USE slepcsys
  USE slepceps
  USE slepcnep
  ! --
  USE solver_context
  IMPLICIT NONE
  CHARACTER(len=128) :: str
  ! --- matrix-free matrix evp (non-linear)
  NEP :: nep
  Mat :: F,J
  TYPE(ApplicationCtx) :: ctx
  TYPE(MatCtx) :: ctxF,ctxJ
  TYPE(MatCtx), POINTER :: ctxF_pt
  PetscInt :: n=128,nev
  KSP :: ksp
  PC :: pc
  PetscMPIInt :: size
  PetscBool :: terse, flg
  PetscErrorCode :: ierr
  ! ---
  EXTERNAL :: FormFunction
  EXTERNAL :: FormJacobian

  ! initialize SLEPc
  CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
  CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,size,ierr)
  IF(size.NE.1) THEN
 CALL MPI_ABORT(PETSC_COMM_WORLD,-1,ierr)
  END IF
  CALL PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-n",n,flg,ierr)
  CALL PetscPrintf(PETSC_COMM_WORLD,"**\n",ierr)
  CALL PetscPrintf(PETSC_COMM_WORLD,"1-D Nonlinear Eigenproblem\n",ierr)
  WRITE(str,*) "n = ",n,"\n"
  CALL PetscPrintf(PETSC_COMM_WORLD,TRIM(str),ierr)
  CALL PetscPrintf(PETSC_COMM_WORLD,"**\n",ierr)
  ctx%h = 1.0d0/n
  ctx%kappa = 1.0;

  ! Create nonlinear eigensolver context
  CALL NEPCreate(PETSC_COMM_WORLD,nep,ierr)

  ! --
  ! Create matrix data structure; set Function evaluation routine
  ctxF%h = ctx%h
  ctxF%kappa = ctx%kappa

  CALL PetscPrintf(PETSC_COMM_WORLD,"Function evaluation routine\n",ierr)

  CALL MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,ctxF,F,ierr)
  CALL MatShellSetOperation(F,MATOP_MULT,MatMult_Fun,ierr)
  CALL MatShellSetOperation(F,MATOP_GET_DIAGONAL,MatGetDiagonal_Fun,ierr)
  CALL MatShellSetOperation(F,MATOP_DESTROY,MatDestroy_Fun,ierr)

  ! Set function matrix data structure and default function evaluation routine
  CALL 

Re: [petsc-users] MatPtAP problem after subsequent call to MatDuplicate

[Samuel Lanthaler]

Thank you for your swift reply, Hong!

Yes, what you said is exactly what I'm doing. No, I don't get an error 
when doing a MatView of matC.


Samuel

On 12/05/2017 04:38 PM, Hong wrote:

Samuel:
You try to do following:
1) Create A;
2) Create P;
3) C = PtAP:
CALL MatPtAP(matA,matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,matC,ierr);
4) MatDestroy(matA,ierr);
5) MatDuplicate(matC,MAT_COPY_VALUES,matA,ierr);
6) MatView(matA,PETSC_VIEWER_STDOUT_WORLD,ierr);

The error occurs at (6). Do you get any error for
MatView(matC)?

C has some special data structures as a matrix product which could be 
lost during

MatDuplicate for matA. It might be a bug in our code. I'll check it.

Hong

Hi there,

I am getting error messages after using MatPtAP to create a new
matrix C = Pt*A*P and then trying to assign A=C. The following is
a minimal example reproducing the problem:

#include "slepc/finclude/slepc.h"
  USE slepcsys
  USE slepceps
  IMPLICIT NONE
  LOGICAL :: cause_error
  ! --- pure PETSc
  PetscErrorCode :: ierr
  PetscScalar :: vals(3,3), val
  Mat :: matA, matP, matC
  PetscInt :: m,idxn(3),idxm(3),idone(1)

  ! initialize SLEPc & Petsc etc.
  CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)

  ! Set up a new matrix
  m = 3
  CALL MatCreate(PETSC_COMM_WORLD,matA,ierr); CHKERRQ(ierr);
  CALL MatSetType(matA,MATMPIAIJ,ierr); CHKERRQ(ierr);
  CALL
MatSetSizes(matA,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr);
CHKERRQ(ierr);
  CALL

MatMPIAIJSetPreallocation(matA,3,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr);
CHKERRQ(ierr);

  ! [ call to MatSetValues to set values of matA]

  ! assemble matrix
  CALL MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);
  CALL MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);

  ! duplicate matrix
  CALL
MatDuplicate(matA,MAT_DO_NOT_COPY_VALUES,matP,ierr);
CHKERRQ(ierr);

 ! [ call to MatSetValues to set values of matP]

  ! assemble matrix
  CALL MatAssemblyBegin(matP,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);
  CALL MatAssemblyEnd(matP,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);

  ! compute C=Pt*A*P
  cause_error = .TRUE. ! set to .TRUE. to cause error,
.FALSE. seems to work fine
  IF(.NOT.cause_error) THEN
 CALL MatDuplicate(matA,MAT_COPY_VALUES,matC,ierr);
CHKERRQ(ierr);
  ELSE
 CALL
MatPtAP(matA,matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,matC,ierr);
CHKERRQ(ierr);
  END IF

  ! destroy matA and duplicate A=C
  CALL MatDestroy(matA,ierr); CHKERRQ(ierr);
  CALL MatDuplicate(matC,MAT_COPY_VALUES,matA,ierr);
CHKERRQ(ierr);

  ! display resulting matrix A
  CALL MatView(matA,PETSC_VIEWER_STDOUT_WORLD,ierr);
CHKERRQ(ierr);

Whether A and P are assigned any values doesn't seem to matter at
all. The error message I'm getting is:

Mat Object: 1 MPI processes
  type: mpiaij
[0]PETSC ERROR:

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation
Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or
-on_error_attach_debugger
[0]PETSC ERROR: or see
http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind

[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and
Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: likely location of problem given in stack below
[0]PETSC ERROR: -  Stack Frames

[0]PETSC ERROR: Note: The EXACT line numbers in the stack are
not available,
[0]PETSC ERROR:   INSTEAD the line number of the start of
the function
[0]PETSC ERROR:   is given.
[0]PETSC ERROR: [0] MatView_MPIAIJ_PtAP line 23
/home/lanthale/Progs/petsc-3.8.2/src/mat/impls/aij/mpi/mpiptap.c
[0]PETSC ERROR: [0] MatView line 949
/home/lanthale/Progs/petsc-3.8.2/src/mat/interface/matrix.c
[0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Signal received
[0]PETSC ERROR: See
http://www.mcs.anl.gov/petsc/documentation/faq.html
 for

[petsc-users] MatPtAP problem after subsequent call to MatDuplicate

[Samuel Lanthaler]

Hi there,

I am getting error messages after using MatPtAP to create a new matrix C 
= Pt*A*P and then trying to assign A=C. The following is a minimal 
example reproducing the problem:


   #include "slepc/finclude/slepc.h"
  USE slepcsys
  USE slepceps
  IMPLICIT NONE
  LOGICAL :: cause_error
  ! --- pure PETSc
  PetscErrorCode :: ierr
  PetscScalar :: vals(3,3), val
  Mat :: matA, matP, matC
  PetscInt :: m,idxn(3),idxm(3),idone(1)

  ! initialize SLEPc & Petsc etc.
  CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)

  ! Set up a new matrix
  m = 3
  CALL MatCreate(PETSC_COMM_WORLD,matA,ierr); CHKERRQ(ierr);
  CALL MatSetType(matA,MATMPIAIJ,ierr); CHKERRQ(ierr);
  CALL MatSetSizes(matA,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr);
   CHKERRQ(ierr);
  CALL
   
MatMPIAIJSetPreallocation(matA,3,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr);
   CHKERRQ(ierr);

  ! [ call to MatSetValues to set values of matA]

  ! assemble matrix
  CALL MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY,ierr);
   CHKERRQ(ierr);
  CALL MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);

  ! duplicate matrix
  CALL MatDuplicate(matA,MAT_DO_NOT_COPY_VALUES,matP,ierr);
   CHKERRQ(ierr);

 ! [ call to MatSetValues to set values of matP]

  ! assemble matrix
  CALL MatAssemblyBegin(matP,MAT_FINAL_ASSEMBLY,ierr);
   CHKERRQ(ierr);
  CALL MatAssemblyEnd(matP,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);

  ! compute C=Pt*A*P
  cause_error = .TRUE. ! set to .TRUE. to cause error, .FALSE.
   seems to work fine
  IF(.NOT.cause_error) THEN
 CALL MatDuplicate(matA,MAT_COPY_VALUES,matC,ierr);
   CHKERRQ(ierr);
  ELSE
 CALL
   MatPtAP(matA,matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,matC,ierr);
   CHKERRQ(ierr);
  END IF

  ! destroy matA and duplicate A=C
  CALL MatDestroy(matA,ierr); CHKERRQ(ierr);
  CALL MatDuplicate(matC,MAT_COPY_VALUES,matA,ierr); CHKERRQ(ierr);

  ! display resulting matrix A
  CALL MatView(matA,PETSC_VIEWER_STDOUT_WORLD,ierr); CHKERRQ(ierr);

Whether A and P are assigned any values doesn't seem to matter at all. 
The error message I'm getting is:


   Mat Object: 1 MPI processes
  type: mpiaij
   [0]PETSC ERROR:
   
   [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation
   Violation, probably memory access out of range
   [0]PETSC ERROR: Try option -start_in_debugger or
   -on_error_attach_debugger
   [0]PETSC ERROR: or see
   http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
   [0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple
   Mac OS X to find memory corruption errors
   [0]PETSC ERROR: likely location of problem given in stack below
   [0]PETSC ERROR: -  Stack Frames
   
   [0]PETSC ERROR: Note: The EXACT line numbers in the stack are not
   available,
   [0]PETSC ERROR:   INSTEAD the line number of the start of the
   function
   [0]PETSC ERROR:   is given.
   [0]PETSC ERROR: [0] MatView_MPIAIJ_PtAP line 23
   /home/lanthale/Progs/petsc-3.8.2/src/mat/impls/aij/mpi/mpiptap.c
   [0]PETSC ERROR: [0] MatView line 949
   /home/lanthale/Progs/petsc-3.8.2/src/mat/interface/matrix.c
   [0]PETSC ERROR: - Error Message
   --
   [0]PETSC ERROR: Signal received
   [0]PETSC ERROR: See
   http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble
   shooting.
   [0]PETSC ERROR: Petsc Release Version 3.8.2, Nov, 09, 2017

Could someone maybe tell me where I'm doing things wrong? The 
documentation for MatPtAP 
 
says that C will be created and that this routine is "currently only 
implemented for pairs of AIJ matrices and classes which inherit from 
AIJ". Does this maybe exclude MPIAIJ matrices?


Thanks,
Samuel

Re: [petsc-users] MatZeroRowsColumns

[Samuel Lanthaler]

Perfect, PETSC_NULL_VEC works in fortran. Thank you very much for your 
help! :-)


Cheers,
Sam

On 11/24/2017 03:18 PM, Matthew Knepley wrote:
On Fri, Nov 24, 2017 at 6:42 AM, Samuel Lanthaler 
<s.lantha...@gmail.com <mailto:s.lantha...@gmail.com>> wrote:


Ah, great. I didn't understand that "optional" arguments are not
optional in the fortran sense of the word. After setting up two
vectors to pass to the routine, the program works. Thank you very
much!

One additional question: In the documentation, it is written that
I can pass PETSC_NULL for these two vectors, if I don't actually
need them. That probably only works in C, but not in fortran. Is
there a corresponding argument I can pass to the routine in
fortran? It seems, that PETSC_NULL_REAL doesn't work; is there
some other PETSC_NULL_XXX that should be passed in this case from
fortran? If not, then I'll just stick to creating vectors and
distroying them afterwards, which is fine for me as well.


You can use PETSC_NULL_VEC

  Thanks,

 Matt

Cheers,
Sam


On 11/24/2017 01:56 AM, Smith, Barry F. wrote:

MatZeroRowsColumns() as two vector arguments you are missing


On Nov 23, 2017, at 2:43 PM, Samuel Lanthaler
<s.lantha...@gmail.com <mailto:s.lantha...@gmail.com>> wrote:

Hi there,
I'm new to PETSc and have been trying to do some basic
manipulation of matrices in fortran, but don't seem to be
able to set a row/column to zero using the
MatZeroRowsColumns command: The following is a small
example program:

! initialize PETSc
   CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)

   ! Set up a new matrix
   m = 3
   CALL MatCreate(PETSC_COMM_WORLD,matA,ierr); CHKERRQ(ierr);
   CALL MatSetType(matA,MATMPIAIJ,ierr); CHKERRQ(ierr);
   CALL
MatSetSizes(matA,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr);
CHKERRQ(ierr);
   CALL

MatMPIAIJSetPreallocation(matA,3,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr);
CHKERRQ(ierr);

   ! set values of matrix
   vals(1,:) = (/1.,2.,3./)
   vals(2,:) = (/4.,5.,6./)
   vals(3,:) = (/7.,8.,9./)
   !
   idxm = (/0,1,2/)
   idxn = (/0,1,2/)
   !
   CALL
MatSetValues(matA,3,idxm,3,idxn,vals,INSERT_VALUES,ierr);
CHKERRQ(ierr);
  ! assemble matrix
   CALL MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);
   CALL MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);
 ! set one row/column to zero, put 6.0d0 on diagonal
   idone(1) = 2
   val = 6.0d0
   CALL MatZeroRowsColumns(matA,1,idone,val,ierr);
CHKERRQ(ierr);

! finalize PETSc
   CALL PetscFinalize(PETSC_NULL_CHARACTER,ierr)

When running the program, I get the following error message:

[0]PETSC ERROR:


[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation
Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or
-on_error_attach_debugger
[0]PETSC ERROR: or see
http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
<http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind>
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux
and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: likely location of problem given in stack
below
[0]PETSC ERROR: -  Stack Frames

[0]PETSC ERROR: Note: The EXACT line numbers in the stack
are not available,
[0]PETSC ERROR:   INSTEAD the line number of the start
of the function
[0]PETSC ERROR:   is given.
[0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Signal received
[0]PETSC ERROR: See
http://www.mcs.anl.gov/petsc/documentation/faq.html
<http://www.mcs.anl.gov/petsc/documentation/faq.html> for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.8.2, Nov, 09, 2017
[0]PETSC ERROR: ./test on a arch-complex-debug named
sam-ThinkPad-T450s by sam Thu Nov 23 23:28:15 2017
[0]PETSC ERROR: Configure options
PETSC_DIR=/home/sam/Progs/petsc-3

Re: [petsc-users] MatZeroRowsColumns

[Samuel Lanthaler]

Ah, great. I didn't understand that "optional" arguments are not 
optional in the fortran sense of the word. After setting up two vectors 
to pass to the routine, the program works. Thank you very much!


One additional question: In the documentation, it is written that I can 
pass PETSC_NULL for these two vectors, if I don't actually need them. 
That probably only works in C, but not in fortran. Is there a 
corresponding argument I can pass to the routine in fortran? It seems, 
that PETSC_NULL_REAL doesn't work; is there some other PETSC_NULL_XXX 
that should be passed in this case from fortran? If not, then I'll just 
stick to creating vectors and distroying them afterwards, which is fine 
for me as well.


Cheers,
Sam


On 11/24/2017 01:56 AM, Smith, Barry F. wrote:

MatZeroRowsColumns() as two vector arguments you are missing



On Nov 23, 2017, at 2:43 PM, Samuel Lanthaler <s.lantha...@gmail.com> wrote:

Hi there,
I'm new to PETSc and have been trying to do some basic manipulation of matrices 
in fortran, but don't seem to be able to set a row/column to zero using the 
MatZeroRowsColumns command: The following is a small example program:

! initialize PETSc
   CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)

   ! Set up a new matrix
   m = 3
   CALL MatCreate(PETSC_COMM_WORLD,matA,ierr); CHKERRQ(ierr);
   CALL MatSetType(matA,MATMPIAIJ,ierr); CHKERRQ(ierr);
   CALL MatSetSizes(matA,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr); CHKERRQ(ierr);
   CALL 
MatMPIAIJSetPreallocation(matA,3,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr); 
CHKERRQ(ierr);

   ! set values of matrix
   vals(1,:) = (/1.,2.,3./)
   vals(2,:) = (/4.,5.,6./)
   vals(3,:) = (/7.,8.,9./)
   !
   idxm = (/0,1,2/)
   idxn = (/0,1,2/)
   !
   CALL MatSetValues(matA,3,idxm,3,idxn,vals,INSERT_VALUES,ierr); CHKERRQ(ierr);
   
   ! assemble matrix

   CALL MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);
   CALL MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);
  
   ! set one row/column to zero, put 6.0d0 on diagonal

   idone(1) = 2
   val = 6.0d0
   CALL MatZeroRowsColumns(matA,1,idone,val,ierr); CHKERRQ(ierr);

! finalize PETSc
   CALL PetscFinalize(PETSC_NULL_CHARACTER,ierr)

When running the program, I get the following error message:

[0]PETSC ERROR: 

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably 
memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see 
http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to 
find memory corruption errors
[0]PETSC ERROR: likely location of problem given in stack below
[0]PETSC ERROR: -  Stack Frames 

[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,
[0]PETSC ERROR:   INSTEAD the line number of the start of the function
[0]PETSC ERROR:   is given.
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: Signal received
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.8.2, Nov, 09, 2017
[0]PETSC ERROR: ./test on a arch-complex-debug named sam-ThinkPad-T450s by sam 
Thu Nov 23 23:28:15 2017
[0]PETSC ERROR: Configure options PETSC_DIR=/home/sam/Progs/petsc-3.8.2 
PETSC_ARCH=arch-complex-debug --with-scalar-type=complex
[0]PETSC ERROR: #1 User provided function() line 0 in  unknown file
--
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 59.
Would someone be so kind as to tell me what I'm doing wrong? Thank you!

Best regards,
Sam

[petsc-users] MatZeroRowsColumns

[Samuel Lanthaler]

Hi there,

I'm new to PETSc and have been trying to do some basic manipulation of 
matrices in fortran, but don't seem to be able to set a row/column to 
zero using the MatZeroRowsColumns command: The following is a small 
example program:|

|

|! initialize PETSc||
||  CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)||
||
||  ! Set up a new matrix||
||  m = 3||
||  CALL MatCreate(PETSC_COMM_WORLD,matA,ierr); CHKERRQ(ierr);||
||  CALL MatSetType(matA,MATMPIAIJ,ierr); CHKERRQ(ierr);||
||  CALL MatSetSizes(matA,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr); 
CHKERRQ(ierr);||
||  CALL 
MatMPIAIJSetPreallocation(matA,3,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr); 
CHKERRQ(ierr);||

||
||  ! set values of matrix||
||  vals(1,:) = (/1.,2.,3./)||
||  vals(2,:) = (/4.,5.,6./)||
||  vals(3,:) = (/7.,8.,9./)||
||  ! ||
||  idxm = (/0,1,2/)||
||  idxn = (/0,1,2/)||
||  !||
||  CALL MatSetValues(matA,3,idxm,3,idxn,vals,INSERT_VALUES,ierr); 
CHKERRQ(ierr);||


||  ! assemble matrix||
||  CALL MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);||
||  CALL MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);||

||  ! set one row/column to zero, put 6.0d0 on diagonal||
||  idone(1) = 2||
||  val = 6.0d0||
||  CALL MatZeroRowsColumns(matA,1,idone,val,ierr); CHKERRQ(ierr);|||

|! finalize PETSc||
||  CALL PetscFinalize(PETSC_NULL_CHARACTER,ierr)||
||
|When running the program, I get the following error message:

|[0]PETSC ERROR: 
||
||[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range||
||[0]PETSC ERROR: Try option -start_in_debugger or 
-on_error_attach_debugger||
||[0]PETSC ERROR: or see 
http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind||
||[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac 
OS X to find memory corruption errors||

||[0]PETSC ERROR: likely location of problem given in stack below||
||[0]PETSC ERROR: -  Stack Frames 
||
||[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not 
available,||
||[0]PETSC ERROR:   INSTEAD the line number of the start of the 
function||

||[0]PETSC ERROR:   is given.||
||[0]PETSC ERROR: - Error Message 
--||

||[0]PETSC ERROR: Signal received||
||[0]PETSC ERROR: See 
http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.||

||[0]PETSC ERROR: Petsc Release Version 3.8.2, Nov, 09, 2017 ||
||[0]PETSC ERROR: ./test on a arch-complex-debug named 
sam-ThinkPad-T450s by sam Thu Nov 23 23:28:15 2017||
||[0]PETSC ERROR: Configure options 
PETSC_DIR=/home/sam/Progs/petsc-3.8.2 PETSC_ARCH=arch-complex-debug 
--with-scalar-type=complex||

||[0]PETSC ERROR: #1 User provided function() line 0 in  unknown file||
||--||
||MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD ||
||with errorcode 59.|

Would someone be so kind as to tell me what I'm doing wrong? Thank you!

Best regards,
Sam