[petsc-users] Fwd: Re: Problem with fortran version of ex29 in ksp
TAY wee-beng
zonexo at gmail.com
Sat Jun 2 04:27:48 CDT 2012
Hi,
I have modified the C version of the ex29 example in ksp to Fortran. I
have attached it. Feel free to include it in the distribution of PETSc.
The DMDAVecGetArrayF90 and DMDAVecRestoreArrayF90 are not working in
linux with ifort, so VecSetValues was used. Worked in VS2008 though.
Btw, I need to solve the Poisson equation multiple times using different
RHS vector.
Which subroutines do I need to call only once?
Can I also add :
call MatSetOption(jac,MAT_STRUCTURALLY_SYMMETRIC,PETSC_TRUE,ierr)
call MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
call MatSetOption(jac,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE,ierr)
call KSPSetOperators(ksp,A_mat,A_mat,SAME_NONZERO_PATTERN,ierr)
Yours sincerely,
TAY wee-beng
On 11/5/2012 4:04 PM, John Mousel wrote:
> Include both header files.
>
> #include <finclude/petscdmda.h>
> #include <finclude/petscdmda.h90>
>
> The h90 header only includes the interface definitions for Fortran.
>
> John
>
> On Fri, May 11, 2012 at 8:51 AM, TAY wee-beng <zonexo at gmail.com
> <mailto:zonexo at gmail.com>> wrote:
>
>
>
> On 11/5/2012 3:30 PM, Matthew Knepley wrote:
>> On Fri, May 11, 2012 at 9:24 AM, TAY wee-beng <zonexo at gmail.com
>> <mailto:zonexo at gmail.com>> wrote:
>>
>> Hi,
>>
>> I have been using the GUI environment to do debugging so I am
>> a bit reluctant to learn Valgrind as its outputs seems a bit
>> daunting. But I guess John is right. I've been spending
>> these few days learning bit by bit.
>>
>> I realised that the error occurs in computerhs, at:
>>
>>
>> I bet this is a beautiful Fortranism. Do you include the F90
>> header file with the interface definition?
>> If not, Fortran just craps out like this. I can't stress enough
>> how much time would be saved by
>> switching languages to something with at least a modicum of error
>> checking.
>
> I initially used:
>
> #include "finclude/petsc.h90"
>
> Compilation and linking was fine in Linux and vs2008.
>
> Now I changed it to what ex22.F was using :
>
> #include <finclude/petscsys.h>
> #include <finclude/petscvec.h>
> #include <finclude/petscmat.h>
> #include <finclude/petscpc.h>
> #include <finclude/petscksp.h>
> #include <finclude/petscdmda.h>
>
> Compiling was ok but linking failed in Linux and VS2008:
>
> undefined reference to `dmdavecgetarrayf90_'
>
> I tried changing #include <finclude/petscdmda.h> to #include
> <finclude/petscdmda.h90> and everything was ok in VS2008 again,
> giving the right answers.
>
> However, in Linux, I got the following error:
>
> [wtay at hpc12:tutorials]$ /opt/openmpi-1.5.3/bin/mpif90 -c -fPIC
> -g  -I/home/wtay/Lib/petsc-3.2-dev_shared_debug/include
> -I/home/wtay/Lib/petsc-3.2-dev_shared_debug/include
> -I/opt/openmpi-1.5.3/include   -o ex29f.o ex29f.F90
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> Â Â Â Â Â Â Â DMDABoundaryType :: pt
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
>        DMDAStencilType :: st
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> Â Â Â Â Â Â Â DMDABoundaryType :: pt
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
>        DMDAStencilType :: st
> --------^
> ex29f.F90(68): error #6404: This name does not have a type, and
> must have an explicit type.  [DMDA_BOUNDARY_NONE]
> call
> DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,i3,i3,PETSC_DECIDE,PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
> -----------------------------------^
> ex29f.F90(68): error #6404: This name does not have a type, and
> must have an explicit type.  [DMDA_STENCIL_STAR]
> call
> DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,i3,i3,PETSC_DECIDE,PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
> -------------------------------------------------------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> Â Â Â Â Â Â Â DMDABoundaryType :: pt
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
>        DMDAStencilType :: st
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> Â Â Â Â Â Â Â DMDABoundaryType :: pt
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
>        DMDAStencilType :: st
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
> Â Â Â Â Â Â Â DMDABoundaryType :: pt
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #5082: Syntax error, found '::' when expecting one of: ( % :
> . = =>
>        DMDAStencilType :: st
> -------------------------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(25):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
> Â Â Â Â Â Â Â DMDABoundaryType :: pt
> --------^
> /home/wtay/Lib/petsc-3.2-dev_shared_debug/include/finclude/ftn-custom/petscdmda.h90(26):
> error #6590: This statement is not permitted as a statement within
> a derived-type-def
>        DMDAStencilType :: st
>
> Is there some errors in petscdmda.h90?
>
>
>
>>
>> Â Â Matt
>> Â
>>
>> *call DMDAVecRestoreArrayF90(da,b,array,ierr)*
>>
>> ==27464== Invalid write of size 8
>> ==27464==Â Â Â at 0x402835: computerhs_ (ex29f.F90:119)
>> ==27464==Â Address 0xfffffffffffffef0 is not stack'd,
>> malloc'd or (recently) free'd
>> ==27464==
>> [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] DM user function line 0 unknownunknown
>> [0]PETSC ERROR: [0] DMComputeFunction line 2085
>> /home/wtay/Codes/petsc-dev/src/dm/interface/dm.c
>> [0]PETSC ERROR: [0] KSPSetUp line 182
>> /home/wtay/Codes/petsc-dev/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: --------------------- Error Message
>> ------------------------------------
>> [0]PETSC ERROR: Signal received!
>>
>> I have checked that "array" 's values are correct. This
>> statement executed without problems in VS2008. If I replace
>> the above with something like:
>>
>> *call VecSet(b,Hy,ierr)*
>>
>> Everything is fine in Linux.
>>
>> Is there something wrong with *DMDAVecRestoreArrayF90*?
>>
>> My code in the area is:
>>
>> call
>> DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
>>
>> call DMDAVecGetArrayF90(da,b,array,ierr)
>>
>> do j = ys,ys+ym-1
>>
>> Â Â Â do i = xs,xs+xm-1
>>
>> Â Â Â Â Â Â array(i,j) =
>> exp(-(i*Hx)*(i*Hx)/nu)*exp(-(j*Hy)*(j*Hy)/nu)*Hx*Hy
>>
>> Â Â Â end do
>>
>> end do
>>
>> call DMDAVecRestoreArrayF90(da,b,array,ierr)
>>
>> call VecAssemblyBegin(b,ierr)
>>
>> call VecAssemblyEnd(b,ierr)
>>
>>
>> Yours sincerely,
>>
>> TAY wee-beng
>>
>>
>> On 8/5/2012 9:41 PM, John Mousel wrote:
>>> TAY wee-bing,
>>>
>>> If you want to be a programmer that writes interesting and
>>> reliable code, you need to be willing to use the tools of
>>> the trade. I can't think of a bigger time-saver than
>>> Valgrind. I would suggest that you learn to use it and use
>>> it a lot. I bet it will lead you to the root of your problem
>>> pretty quickly.
>>>
>>> John
>>>
>>> On Tue, May 8, 2012 at 2:17 PM, TAY wee-beng
>>> <zonexo at gmail.com <mailto:zonexo at gmail.com>> wrote:
>>>
>>> Hi,
>>>
>>> I compiled and run my code under visual studio 2008 with
>>> intel fortran. Everything works ok.
>>>
>>> However, when I tried to run the code in linux, I got
>>> the error as below. The error happens when
>>> KSPSetUp(ksp,ierr) is called.
>>>
>>> However, I am not able to print VecView or MatView to
>>> view if there's any errors. Is there any recommendation
>>> for debugging? I hope I do not need to valgrind if possible.
>>>
>>> [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] DM user function line 0 unknownunknown
>>> [0]PETSC ERROR: [0] DMComputeFunction line 2085
>>> /home/wtay/Codes/petsc-dev/src/dm/interface/dm.c
>>> [0]PETSC ERROR: [0] KSPSetUp line 182
>>> /home/wtay/Codes/petsc-dev/src/ksp/ksp/interface/itfunc.c
>>> [0]PETSC ERROR: --------------------- Error Message
>>> ------------------------------------
>>> [0]PETSC ERROR: Signal received!
>>> [0]PETSC ERROR:
>>> ------------------------------------------------------------------------
>>> [0]PETSC ERROR: Petsc Development HG revision:
>>> 7ecdd63ec420b1659b960e65d96e822c5ac1a968Â HG Date: Mon
>>> May 07 21:42:26 2012 -0500
>>> [0]PETSC ERROR: See docs/changes/index.html for recent
>>> updates.
>>> [0]PETSC ERROR: See docs/faq.html for hints about
>>> trouble shooting.
>>> [0]PETSC ERROR: See docs/index.html for manual pages.
>>> [0]PETSC ERROR:
>>> ------------------------------------------------------------------------
>>> [0]PETSC ERROR: ./ex29f on a petsc-3.2 named hpc12 by
>>> wtay Tue May 8 20:45:42 2012
>>> [0]PETSC ERROR: Libraries linked from
>>> /home/wtay/Lib/petsc-3.2-dev_shared_debug/lib
>>> [0]PETSC ERROR: Configure run at Tue May 8 10:47:59 2012
>>> [0]PETSC ERROR: Configure options
>>> --with-mpi-dir=/opt/openmpi-1.5.3/
>>> --with-blas-lapack-dir=/opt/intelcpro-11.1.059/mkl/lib/em64t/
>>> --with-debugging=1 --download-hypre=1
>>> --prefix=/home/wtay/Lib/petsc-3.2-dev_shared_debug
>>> --known-mpi-shared=1 --with-shared-libraries
>>> [0]PETSC ERROR:
>>> ------------------------------------------------------------------------
>>> [0]PETSC ERROR: User provided function() line 0 in
>>> unknown directory unknown file
>>> --------------------------------------------------------------------------
>>> MPI_ABORT was invoked on rank 0 in communicator
>>> MPI_COMM_WORLD
>>> with errorcode 59.
>>>
>>> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI
>>> processes.
>>> You may or may not see output from other processes,
>>> depending on
>>> exactly when Open MPI kills them.
>>>
>>> Yours sincerely,
>>>
>>> TAY wee-beng
>>>
>>>
>>> On 5/5/2012 1:43 AM, Matthew Knepley wrote:
>>>> On Fri, May 4, 2012 at 5:42 PM, TAY wee-beng
>>>> <zonexo at gmail.com <mailto:zonexo at gmail.com>> wrote:
>>>>
>>>> Hi,
>>>>
>>>> I wonder if you people are interested to include my
>>>> ex29 fortran version in the petsc examples, which
>>>> can help people who are using fortran.
>>>>
>>>>
>>>> Yes, we will definitely include it. Please send the
>>>> source, and a representative output with run options.
>>>>
>>>> Â Thanks,
>>>>
>>>> Â Â Â Matt
>>>> Â
>>>>
>>>> Thanks.
>>>>
>>>> Yours sincerely,
>>>>
>>>> TAY wee-beng
>>>>
>>>>
>>>> On 4/5/2012 9:28 PM, Matthew Knepley wrote:
>>>>> On Fri, May 4, 2012 at 3:24 PM, TAY wee-beng
>>>>> <zonexo at gmail.com <mailto:zonexo at gmail.com>> wrote:
>>>>>
>>>>>
>>>>> On 4/5/2012 9:16 PM, Barry Smith wrote:
>>>>>
>>>>> Â Â Do an hg pull and then run make in
>>>>> src/mat/interface/ftn-custom/
>>>>>
>>>>> Â Â Then it should link.
>>>>>
>>>>> Â Â Barry
>>>>>
>>>>> Â Â There was a E missing from the all
>>>>> caps name of the function.
>>>>>
>>>>> After hg pull, I did:
>>>>>
>>>>> cd src//mat/interface/ftn-custom/
>>>>>
>>>>> User at User-PC
>>>>> /cygdrive/c/temp/petsc-dev/src/mat/interface/ftn-custom
>>>>> $ make
>>>>> make[1]: Warning: File
>>>>> `/cygdrive/c/temp/petsc-dev/petsc-3.2-dev_win32_vs2008/lib/libpetsc.lib(zmatregf.o)'
>>>>> has modification time 787 s in the future
>>>>> make[1]: Nothing to be done for `libc'.
>>>>> make[1]: warning: Â Clock skew detected.
>>>>> Â Your build may be incomplete.
>>>>>
>>>>> But it still can't work.
>>>>>
>>>>>
>>>>> Something is messed up with the clock on this machine.
>>>>>
>>>>> HOWEVER, development requires certain basic skills
>>>>> in order to debug your work. We
>>>>> cannot be the ones debugging your code. Now
>>>>>
>>>>> Â nm $PETSC_ARCH/lib/libpetsc.a | grep -i
>>>>> MatNullSpaceRemove
>>>>>
>>>>> Â will look for the symbol.
>>>>>
>>>>> Â Â Matt
>>>>> Â
>>>>>
>>>>>
>>>>>
>>>>> On May 4, 2012, at 2:11 PM, Matthew
>>>>> Knepley wrote:
>>>>>
>>>>> On Fri, May 4, 2012 at 3:01 PM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>> Â wrote:
>>>>>
>>>>> On 4/5/2012 5:17 PM, Matthew Knepley
>>>>> wrote:
>>>>>
>>>>> On Fri, May 4, 2012 at 11:05 AM,
>>>>> TAY wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>> Â wrote:
>>>>>
>>>>> On 4/5/2012 3:05 PM, Matthew
>>>>> Knepley wrote:
>>>>>
>>>>> On Fri, May 4, 2012 at 8:59
>>>>> AM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>>
>>>>> Â wrote:
>>>>>
>>>>> Hi,
>>>>>
>>>>> Is there anything else I can
>>>>> try to get it working right?
>>>>>
>>>>> The MatGetNullSpaceRemove() is
>>>>> missing.
>>>>>
>>>>> Where should I add
>>>>> MatGetNullSpaceRemove and what are
>>>>> its syntax? I googled but there's
>>>>> no results.
>>>>>
>>>>> Fixed in p;etsc-dev:
>>>>>
>>>>> Â
>>>>> http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MatNullSpaceRemove.html
>>>>>
>>>>> I just compiled the updated petsc-dev
>>>>> but I got the same error msg when I use:
>>>>>
>>>>> call
>>>>> MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)
>>>>>
>>>>> error LNK2019: unresolved external
>>>>> symbol MATNULLSPACEREMOVE referenced
>>>>> in function COMPUTERHS
>>>>> 1>c:\obj_tmp\ex29f\Debug\ex29f.exe :
>>>>> fatal error LNK1120: 1 unresolved
>>>>> externals
>>>>>
>>>>> That function is in:
>>>>>
>>>>> Â src/mat/interface/ftn-custom/zmatrixf.c
>>>>>
>>>>> Did that compile? Can you see the
>>>>> symbol in libpetsc.a?
>>>>>
>>>>> Â Â Matt
>>>>>
>>>>> Â Â Â Matt
>>>>>
>>>>> Thanks.
>>>>>
>>>>> Â Â Matt
>>>>>
>>>>> Thanksï¼
>>>>>
>>>>>
>>>>> On 2/5/2012 10:11 PM, Matthew
>>>>> Knepley wrote:
>>>>>
>>>>> On Wed, May 2, 2012 at
>>>>> 1:55 PM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>>
>>>>> Â wrote:
>>>>> Hi,
>>>>>
>>>>> I did a MatView and
>>>>> VecView on both C and
>>>>> Fortran, right after Mat
>>>>> and Vec assembly. I have
>>>>> attached the printout
>>>>> below. They are exactly
>>>>> the same, but yet the
>>>>> result is different in
>>>>> Neumann condition.
>>>>> However, the dirichlet
>>>>> condition gives the
>>>>> correct ans. Is there
>>>>> anything else that could
>>>>> be wrong even if the Mat
>>>>> and Vec are the same?
>>>>>
>>>>> Did you set the null space
>>>>> for the matrix when you
>>>>> have Neumann conditions?
>>>>>
>>>>> Yes, for the matrix, I set as:
>>>>>
>>>>> call
>>>>> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
>>>>>
>>>>> Â Â call
>>>>> MatSetNullSpace(jac,nullspace,ierr)
>>>>>
>>>>> Â Â call
>>>>> MatNullSpaceDestroy(nullspace,ierr)
>>>>>
>>>>> for the Vec,
>>>>>
>>>>> call
>>>>> MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
>>>>>
>>>>> Â Â !call
>>>>> MatNullSpaceRemove(nullspace,b,PETSC_NULL,ierr)
>>>>>
>>>>> Â Â call
>>>>> MatNullSpaceDestroy(nullspace,ierr)
>>>>>
>>>>> MatNullSpaceRemove was comment
>>>>> out because there's error
>>>>> during linking
>>>>>
>>>>>
>>>>> Â Â Matt
>>>>>
>>>>> Thanks!
>>>>>
>>>>> Fortran:
>>>>>
>>>>> Matrix Object: 1 MPI processes
>>>>> Â type: seqaij
>>>>> row 0: (0, 2) Â (1, -1)
>>>>> Â (3, -1)
>>>>> row 1: (0, -1) Â (1, 3)
>>>>> Â (2, -1) Â (4, -1)
>>>>> row 2: (1, -1) Â (2, 2)
>>>>> Â (5, -1)
>>>>> row 3: (0, -1) Â (3, 3)
>>>>> Â (4, -1) Â (6, -1)
>>>>> row 4: (1, -1) Â (3, -1)
>>>>> Â (4, 4) Â (5, -1) Â (7, -1)
>>>>> row 5: (2, -1) Â (4, -1)
>>>>> Â (5, 3) Â (8, -1)
>>>>> row 6: (3, -1) Â (6, 2)
>>>>> Â (7, -1)
>>>>> row 7: (4, -1) Â (6, -1)
>>>>> Â (7, 3) Â (8, -1)
>>>>> row 8: (5, -1) Â (7, -1)
>>>>> Â (8, 2)
>>>>> Vector
>>>>> Object:Vec_0000000084000000_0
>>>>> 1 MPI processes
>>>>> Â type: mpi
>>>>> Process [0]
>>>>> 0.25
>>>>> 0.0205213
>>>>> 1.135e-005
>>>>> 0.0205213
>>>>> 0.00168449
>>>>> 9.31663e-007
>>>>> 1.135e-005
>>>>> 9.31663e-007
>>>>> 5.15289e-010
>>>>> Vector
>>>>> Object:Vec_0000000084000000_1
>>>>> 1 MPI processes
>>>>> Â type: mpi
>>>>> Process [0]
>>>>> 0.14924
>>>>> 0.0242397
>>>>> -0.0260347
>>>>> 0.0242397
>>>>> -0.0256192
>>>>> -0.0400102
>>>>> -0.0260347
>>>>> -0.0400102
>>>>> -0.0400102
>>>>> Press any key to continue
>>>>> . . .
>>>>>
>>>>> C:
>>>>>
>>>>> Matrix Object: 1 MPI processes
>>>>> Â type: seqaij
>>>>> row 0: (0, 2) Â (1, -1)
>>>>> Â (3, -1)
>>>>> row 1: (0, -1) Â (1, 3)
>>>>> Â (2, -1) Â (4, -1)
>>>>> row 2: (1, -1) Â (2, 2)
>>>>> Â (5, -1)
>>>>> row 3: (0, -1) Â (3, 3)
>>>>> Â (4, -1) Â (6, -1)
>>>>> row 4: (1, -1) Â (3, -1)
>>>>> Â (4, 4) Â (5, -1) Â (7, -1)
>>>>> row 5: (2, -1) Â (4, -1)
>>>>> Â (5, 3) Â (8, -1)
>>>>> row 6: (3, -1) Â (6, 2)
>>>>> Â (7, -1)
>>>>> row 7: (4, -1) Â (6, -1)
>>>>> Â (7, 3) Â (8, -1)
>>>>> row 8: (5, -1) Â (7, -1)
>>>>> Â (8, 2)
>>>>> Vector
>>>>> Object:Vec_0x1d3b000_0 1
>>>>> MPI processes
>>>>> Â type: mpi
>>>>> Process [0]
>>>>> 0.25
>>>>> 0.0205212
>>>>> 1.135e-05
>>>>> 0.0205212
>>>>> 0.00168449
>>>>> 9.31663e-07
>>>>> 1.135e-05
>>>>> 9.31663e-07
>>>>> 5.15288e-10
>>>>> Vector
>>>>> Object:Vec_0x1d3b000_1 1
>>>>> MPI processes
>>>>> Â type: mpi
>>>>> Process [0]
>>>>> 0.139311
>>>>> 0.0305751
>>>>> -0.0220633
>>>>> 0.0305751
>>>>> -0.0135158
>>>>> -0.042185
>>>>> -0.0220633
>>>>> -0.042185
>>>>> -0.058449
>>>>>
>>>>>
>>>>>
>>>>> Yours sincerely,
>>>>>
>>>>> TAY wee-beng
>>>>>
>>>>>
>>>>> On 1/5/2012 11:54 PM,
>>>>> Matthew Knepley wrote:
>>>>>
>>>>> On Tue, May 1, 2012 at
>>>>> 5:48 PM, TAY
>>>>> wee-beng<zonexo at gmail.com
>>>>> <mailto:zonexo at gmail.com>>
>>>>> Â wrote:
>>>>> Hi,
>>>>>
>>>>> Do you mean my method
>>>>> is wrong?
>>>>>
>>>>> I am following the
>>>>> template of ex22f,
>>>>>
>>>>> where the variables
>>>>> are declared as :
>>>>>
>>>>> PetscScalar  v(5)
>>>>>
>>>>> MatStencil Â
>>>>> row(4),col(4,5)
>>>>>
>>>>> Hence,
>>>>>
>>>>> for the neumann BC
>>>>>
>>>>> num = 1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â if (j/=0) then
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â v(num) =
>>>>> -rho*HxdHy
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_i,num) =
>>>>> i
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_j,num) =
>>>>> j-1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â num = num + 1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â end if
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â if (i/=0) then
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â v(num) =
>>>>> -rho*HydHx
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_i,num) =
>>>>> i-1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_j,num) =
>>>>> j
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â num = num + 1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â end if
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â if (i/=mx-1) then
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â v(num) =
>>>>> -rho*HydHx
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_i,num) =
>>>>> i+1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_j,num) =
>>>>> j
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â num = num + 1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â end if
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â if (j/=my-1) then
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â v(num) =
>>>>> -rho*HxdHy
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_i,num) =
>>>>> i
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â
>>>>> Â col(MatStencil_j,num) =
>>>>> j+1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â Â Â num = num + 1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â end if
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â v(num) =
>>>>> ((num-1)/2.0)*rho*(HxdHy
>>>>> + HydHx)
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â print *, v
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â col(MatStencil_i,num) =
>>>>> i
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â col(MatStencil_j,num) =
>>>>> j
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â !num = num + 1
>>>>>
>>>>> Â Â Â Â Â Â Â
>>>>> Â call
>>>>> MatSetValuesStencil(jac,i1,row,num,col,v,INSERT_VALUES,ierr)
>>>>>
>>>>> I do not get any more
>>>>> out of range error.
>>>>> However,my ans is
>>>>> still different from
>>>>> that of ex29 in C.
>>>>>
>>>>> This is very simple.
>>>>> You have an error in
>>>>> your code. Checking it
>>>>> is very simple: run
>>>>> the code and
>>>>> break in
>>>>> MatSetValues(). Make
>>>>> sure ex29 makes calls
>>>>> with exactly the same
>>>>> indices as your ex29f.
>>>>>
>>>>> Â Â Matt
>>>>>
>>>>> Yours sincerely,
>>>>>
>>>>> TAY wee-beng
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> 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
>>>
>>>
>>
>>
>>
>> --
>> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120602/58430b8e/attachment-0001.html>
-------------- next part --------------
! Concepts: KSP^solving a system of linear equations
! Concepts: KSP^Laplacian, 2d
! Processors: n
!Added by Wee Beng Tay (W.B.Tay at tudelft.nl).
!conversion from ex29 C version to Fortran
!Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
! -div \rho grad u = f, 0 < x,y < 1,
!with forcing function
! f = e^{-x^2/\nu} e^{-y^2/\nu}
!with Dirichlet boundary conditions
! u = f(x,y) for x = 0, x = 1, y = 0, y = 1
!or pure Neumman boundary conditions
!This uses multigrid to solve the linear system
program ex29f
implicit none
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Include files
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! petscsys.h - base PETSc routines petscvec.h - vectors
! petscmat.h - matrices
! petscksp.h - Krylov subspace methods petscpc.h - preconditioners
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscpc.h>
#include <finclude/petscksp.h>
#include <finclude/petscdmda.h>
#include <finclude/petscdmda.h90>
!#include "finclude/petsc.h90"
PetscErrorCode ierr
DM da
KSP ksp
Vec x
external ComputeRHS,ComputeMatrix
PetscInt i1,i3,bc_cond !>bc_cond = boundary condition = 1/2 = dirichlet/neumann
!>have to change bc_cond at 3 locations
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
bc_cond = 2
i3 = -3
i1 = 1
call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
call DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,i3,i3,PETSC_DECIDE,PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
call DMSetFunction(da,ComputeRHS,ierr)
call DMSetJacobian(da,ComputeMatrix,ierr)
call KSPSetDM(ksp,da,ierr)
call KSPSetFromOptions(ksp,ierr)
call KSPSetUp(ksp,ierr)
call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
call KSPGetSolution(ksp,x,ierr)
call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)
call KSPDestroy(ksp,ierr)
call DMDestroy(da,ierr)
call PetscFinalize(ierr)
end program ex29f
subroutine ComputeRHS(da,x,b,ierr)
implicit none
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscpc.h>
#include <finclude/petscksp.h>
#include <finclude/petscdmda.h90>
#include <finclude/petscdmda.h>
PetscErrorCode ierr
PetscInt i,j,ij,mx,my,xm,ym,xs,ys,bc_cond,status
PetscScalar h,nu,rho
PetscScalar Hx,Hy
!PetscScalar,pointer :: array(:,:) !Problem with using DMDAVecGetArrayF90 and DMDAVecRestoreArrayF90 below in
Vec x,b !linux with ifort, so VecSetValues was used. Worked in VS2008 though
DM da
MatNullSpace nullspace !>neumann BC
integer, allocatable:: ij_array(:)
real(8), allocatable :: array(:)
nu = 0.1
bc_cond = 2
call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
Hx = 1.d0 / (mx-1)
Hy = 1.d0 / (my-1)
call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
allocate (ij_array(0:(ym-ys)*(xm-xs)-1), STAT=status)
if (status/=0) STOP "Cannot allocate memory"
allocate (array(0:(ym-ys)*(xm-xs)-1), STAT=status)
if (status/=0) STOP "Cannot allocate memory"
!call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,(ym-ys)*(xm-xs),b,ierr)
!call DMDAVecGetArrayF90(da,b,array,ierr)
!do j = ys,ys+ym-1
! do i = xs,xs+xm-1
! array(i,j) = exp(-(i*Hx)*(i*Hx)/nu)*exp(-(j*Hy)*(j*Hy)/nu)*Hx*Hy
! end do
!end do
!call DMDAVecRestoreArrayF90(da,b,array,ierr)
do j = ys,ys+ym-1
do i = xs,xs+xm-1
ij = (j)*(xm-xs) + i
ij_array(ij) = ij
array(ij) = exp(-(i*Hx)*(i*Hx)/nu)*exp(-(j*Hy)*(j*Hy)/nu)*Hx*Hy
end do
end do
call VecSetValues(b,(ym-ys)*(xm-xs),ij_array,array,INSERT_VALUES,ierr)
call VecAssemblyBegin(b,ierr)
call VecAssemblyEnd(b,ierr)
!call VecView(b,PETSC_VIEWER_STDOUT_WORLD,ierr)
if (bc_cond == 2) then
call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
call MatNullSpaceRemove(nullspace,b,PETSC_NULL_OBJECT,ierr)
call MatNullSpaceDestroy(nullspace,ierr)
end if
end subroutine ComputeRHS
subroutine ComputeMatrix(da,x,JJ,jac,str,ierr)
implicit none
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscdmda.h>
#include <finclude/petscdmda.h90>
Mat jac,JJ
PetscErrorCode ierr
DM da
PetscInt i,j,k,mx,my,xm,num
PetscInt ym,xs,ys,i1,i5,bc_cond
PetscScalar v(5),Hx,Hy,rho,centerRho
PetscScalar HydHx
PetscScalar HxdHy
MatStencil row(4),col(4,5)
Vec x
MatStructure str
MatNullSpace nullspace !>neumann BC
bc_cond = 2
rho = 1.0
i1 = 1
i5 = 5
centerRho = rho
call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
Hx = 1.d0 / (mx-1)
Hy = 1.d0 / (my-1)
HxdHy = Hx/Hy
HydHx = Hy/Hx
call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
do j=ys,ys+ym-1
do i=xs,xs+xm-1
row(MatStencil_i) = i
row(MatStencil_j) = j
call ComputeRho(i,j,mx,my,centerRho,rho)
if (i==0 .or. j==0 .or. i==mx-1 .or. j==my-1) then
if (bc_cond == 1) then
v(1) = 2.0*rho*(HxdHy + HydHx)
call MatSetValuesStencil(jac,i1,row,i1,row,v,INSERT_VALUES,ierr)
else if (bc_cond == 2) then
num = 1
if (j/=0) then
v(num) = -rho*HxdHy
col(MatStencil_i,num) = i
col(MatStencil_j,num) = j-1
num = num + 1
end if
if (i/=0) then
v(num) = -rho*HydHx
col(MatStencil_i,num) = i-1
col(MatStencil_j,num) = j
num = num + 1
end if
if (i/=mx-1) then
v(num) = -rho*HydHx
col(MatStencil_i,num) = i+1
col(MatStencil_j,num) = j
num = num + 1
end if
if (j/=my-1) then
v(num) = -rho*HxdHy
col(MatStencil_i,num) = i
col(MatStencil_j,num) = j+1
num = num + 1
end if
v(num) = ((num-1)/2.0)*rho*(HxdHy + HydHx)
col(MatStencil_i,num) = i
col(MatStencil_j,num) = j
!num = num + 1
call MatSetValuesStencil(jac,i1,row,num,col,v,INSERT_VALUES,ierr)
end if
else
v(1) = -rho*HxdHy
col(MatStencil_i,1) = i
col(MatStencil_j,1) = j-1
v(2) = -rho*HydHx
col(MatStencil_i,2) = i-1
col(MatStencil_j,2) = j
v(3) = 2.0*rho*(HxdHy + HydHx)
col(MatStencil_i,3) = i
col(MatStencil_j,3) = j
v(4) = -rho*HydHx
col(MatStencil_i,4) = i+1
col(MatStencil_j,4) = j
v(5) = -rho*HxdHy
col(MatStencil_i,5) = i
col(MatStencil_j,5) = j+1
call MatSetValuesStencil(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)
end if
end do
end do
call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
!call MatView(jac,PETSC_VIEWER_STDOUT_WORLD,ierr)
if (bc_cond == 2) then
call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL_OBJECT,nullspace,ierr)
call MatSetNullSpace(jac,nullspace,ierr)
call MatNullSpaceDestroy(nullspace,ierr)
end if
end subroutine ComputeMatrix
subroutine ComputeRho(i,j,mx,my,centerRho,rho)
PetscInt i,j,mx,my
PetscScalar rho,centerRho
if ((i > mx/3.0) .and. (i < 2.0*mx/3.0) .and. (j > my/3.0) .and. (j < 2.0*my/3.0)) then
rho = centerRho
else
rho = 1.0
end if
end subroutine ComputeRho
More information about the petsc-users
mailing list