<div dir="ltr">It does look like memory is getting corrupted and it just fails first in MatGetSizes. This happens after many calls the same solver routine that does not have any side effects as I recall. Does fail with a different PC like Jacobi? As I recall you checked for memory leaks. You mentioned that it failed at different times when you reordered the Destroy calls. Is this repeatable? (ie, does it fail on the same solve number every time with a given configuration. Could you try commenting out half the Destroys and see if that fixes it (if you don't run out of memory). Then try the other half of the Destroys and if you get a fix go for bisection search.<div>
Mark</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Thu, Nov 21, 2013 at 12:17 PM, Satish Balay <span dir="ltr"><<a href="mailto:balay@mcs.anl.gov" target="_blank">balay@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">matsetsizes() is used practically everywhere - so I can't believe<br>
there is a bug there. [the fortran interface code autogenerated by<br>
bfort].<br>
<br>
Do you see the error with the 'first' call to this routine - or some<br>
subsequent loop?<br>
<br>
Perhaps you can try using valgrind equivalent on windows to see if<br>
there is memory corruption somewhere else.<br>
<br>
<a href="http://stackoverflow.com/questions/413477/is-there-a-good-valgrind-substitute-for-windows" target="_blank">http://stackoverflow.com/questions/413477/is-there-a-good-valgrind-substitute-for-windows</a><br>
suggests using <a href="http://code.google.com/p/drmemory/" target="_blank">http://code.google.com/p/drmemory/</a><br>
<br>
Satish<br>
<div><div class="h5"><br>
On Thu, 21 Nov 2013, Einar Sørheim wrote:<br>
<br>
> I am on the windows platform so unfortunatley no valgrind, however have<br>
> been able to track down the problem in the debugger and it seems like a<br>
> matrixpointer gets corrupted when using the fortran-c interface in petsc.<br>
> At some point in the calculation the matrix pointer in MatSetSize gets<br>
> corrupted when called from fortran, propably in this cast:<br>
><br>
> void PETSC_STDCALL matsetsizes_(Mat A,PetscInt *m,PetscInt *n,PetscInt<br>
> *M,PetscInt *N, int *__ierr ){<br>
> *__ierr = MatSetSizes(<br>
> (Mat)PetscToPointer((A) ),*m,*n,*M,*N);<br>
> }<br>
> when entering the MatSetSize in the debugger the A-matrix pointer is<br>
> corrupt and I get an access violation on code line:<br>
> PetscValidHeaderSpecific(A,MAT_CLASSID,1);<br>
> Any ideas what may be the cause of this?<br>
> Regards<br>
> Einar<br>
><br>
><br>
> 2013/10/24 Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>><br>
><br>
> ><br>
> > So apparently no memory leak.<br>
> ><br>
> > If you are still getting crashes you will need to run the case that<br>
> > crashes 1) with valgrind if you haven’t or 2) using<br>
> > -on_error_attach_debugger or -start_in_debugger to track down the cause of<br>
> > the crash.<br>
> ><br>
> > Barry<br>
> ><br>
> ><br>
> ><br>
> > On Oct 24, 2013, at 2:09 AM, Einar Sørheim <<a href="mailto:einar.sorheim@gmail.com">einar.sorheim@gmail.com</a>><br>
> > wrote:<br>
> ><br>
> > > Well i I did a rerun without -malloc_log, this is what I got (I have<br>
> > also attached the log-sumary file):<br>
> > ><br>
> > > [0] PetscFinalize(): PetscFinalize() called<br>
> > > [0] PetscCommDuplicate(): Duplicating a communicator 1140850688<br>
> > -2080374784 max tags = <a href="tel:2147483647" value="+12147483647">2147483647</a><br>
> > > Summary of Memory Usage in PETSc<br>
> > > [0]Current space PetscMalloc()ed 45936, max space PetscMalloced()<br>
> > 1.13728e+009<br>
> > > [0]OS cannot compute process memory<br>
> > > [0] PetscCommDuplicate(): Using internal PETSc communicator 1140850688<br>
> > -2080374784<br>
> > > [0] PetscGetHostName(): Rejecting domainname, likely is NIS CMP3.<br>
> > > [0] PetscGetHostName(): Rejecting domainname, likely is NIS CMP3.<br>
> > > [0] PetscFOpen(): Opening file Log.0<br>
> > > [0] Petsc_DelViewer(): Removing viewer data attribute in an MPI_Comm<br>
> > -2080374784<br>
> > > [0] Petsc_DelComm_Inner(): Removing reference to PETSc communicator<br>
> > embedded in a user MPI_Comm -2080374784<br>
> > > [0] Petsc_DelComm_Outer(): User MPI_Comm 1140850688 is being freed after<br>
> > removing reference from inner PETSc comm to this outer comm<br>
> > > [0] PetscCommDestroy(): Deleting PETSc MPI_Comm -2080374784<br>
> > > [0] Petsc_DelViewer(): Removing viewer data attribute in an MPI_Comm<br>
> > -2080374784<br>
> > > [0] Petsc_DelThreadComm(): Deleting thread communicator data in an<br>
> > MPI_Comm -2080374784<br>
> > > [0] Petsc_DelCounter(): Deleting counter data in an MPI_Comm -2080374784<br>
> > ><br>
> > ><br>
> > ><br>
> > > 2013/10/23 Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>><br>
> > ><br>
> > > Don’t run with the -malloc_log Mark did not suggest using that option.<br>
> > ><br>
> > ><br>
> > > On Oct 23, 2013, at 8:48 AM, Einar Sørheim <<a href="mailto:einar.sorheim@gmail.com">einar.sorheim@gmail.com</a>><br>
> > wrote:<br>
> > ><br>
> > > > I tried running with the following options :<br>
> > > > -info<br>
> > > > -malloc<br>
> > > > -malloc_log<br>
> > > > -malloc_dump<br>
> > > > -malloc_debug<br>
> > > > -memory_info<br>
> > > > -log_summary log_summary.txt<br>
> > > > -log<br>
> > > > but didn't get a lot wiser however in my case would it be more<br>
> > efficient to keep the ksp solver object, if so should I have one object eg.<br>
> > for the thermal problem and another one for the mechanical?<br>
> > > > The output from the end of the run:<br>
> > > > Total wall clock time actually solving: 856.231474<br>
> > > > [0] PetscFinalize(): PetscFinalize() called<br>
> > > > [0] PetscCommDuplicate(): Duplicating a communicator 1140850688<br>
> > -2080374784 max tags = <a href="tel:2147483647" value="+12147483647">2147483647</a><br>
> > > > Summary of Memory Usage in PETSc<br>
> > > > [0]Current space PetscMalloc()ed 45936, max space PetscMalloced()<br>
> > 1.13728e+009<br>
> > > > [0]OS cannot compute process memory<br>
> > > > [0] PetscCommDuplicate(): Using internal PETSc communicator 1140850688<br>
> > -2080374784<br>
> > > > [0] PetscGetHostName(): Rejecting domainname, likely is NIS CMP3.<br>
> > > > [0] PetscGetHostName(): Rejecting domainname, likely is NIS CMP3.<br>
> > > > [0] PetscFOpen(): Opening file Log.0<br>
> > > > [0] Petsc_DelViewer(): Removing viewer data attribute in an MPI_Comm<br>
> > -2080374784<br>
> > > > [0] Petsc_DelComm_Inner(): Removing reference to PETSc communicator<br>
> > embedded in a user MPI_Comm -2080374784<br>
> > > > [0] Petsc_DelComm_Outer(): User MPI_Comm 1140850688 is being freed<br>
> > after removing reference from inner PETSc comm to this outer comm<br>
> > > > [0] PetscCommDestroy(): Deleting PETSc MPI_Comm -2080374784<br>
> > > > [0] Petsc_DelViewer(): Removing viewer data attribute in an MPI_Comm<br>
> > -2080374784<br>
> > > > [0] Petsc_DelThreadComm(): Deleting thread communicator data in an<br>
> > MPI_Comm -2080374784<br>
> > > > [0] Petsc_DelCounter(): Deleting counter data in an MPI_Comm<br>
> > -2080374784<br>
> > > > [0]PETSC ERROR: --------------------- Error Message<br>
> > ------------------------------------<br>
> > > > [0]PETSC ERROR: Object is in wrong state!<br>
> > > > [0]PETSC ERROR: PetscMallocDumpLog() called without call to<br>
> > PetscMallocSetDumpLog() this is often due to<br>
> > > > setting the option -malloc_log AFTER<br>
> > PetscInitialize() with PetscOptionsInsert() or PetscOptionsInsertFile()!<br>
> > > > [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > > > [0]PETSC ERROR: Petsc Release Version 3.4.3, Oct, 15, 2013<br>
> > > > [0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>
> > > > [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>
> > > > [0]PETSC ERROR: See docs/index.html for manual pages.<br>
> > > > [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > > > [0]PETSC ERROR: D:\Programming\gits\Alsim\SharedModules -<br>
> > Copy\Source\Release\Alsim_nypetsc.exe on a arch-mswin-c-debug named CMP3 by<br>
> > admeinar Wed Oct 23 11:54:03 2013<br>
> > > > [0]PETSC ERROR: Libraries linked from<br>
> > /cygdrive/d/Programming/petsc-3.4.3/arch-mswin-c-debug/lib<br>
> > > > [0]PETSC ERROR: Configure run at Thu Oct 17 08:23:18 2013<br>
> > > > [0]PETSC ERROR: Configure options --with-debugging=1 --with-openmp=yes<br>
> > --with-pthread=no --with-cc="win32fe icl" --with-fc="win32fe ifort"<br>
> > --with-cxx="win32fe icl" --with-blas-lapack-dir="C:/Program Files<br>
> > (x86)/Intel/Composer XE 2013/mkl/lib/intel64" --download-superlu<br>
> > --with-sowing=0 --with-c2html=0 --with-mpi-dir="C:/Program Files<br>
> > (x86)/Intel/MPI/<a href="http://4.1.0.028/em64t" target="_blank">4.1.0.028/em64t</a>" --useThreads=0 --useThreads=0<br>
> > > > [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > > > [0]PETSC ERROR: PetscMallocDumpLog() line 652 in<br>
> > src/sys/memory/D:\PROGRA~1\PETSC-~1.3\src\sys\memory\mtr.c<br>
> > > > [0]PETSC ERROR: PetscFinalize() line 1185 in<br>
> > src/sys/objects/D:\PROGRA~1\PETSC-~1.3\src\sys\objects\pinit.c<br>
> > > ><br>
> > > ><br>
> > > > 2013/10/19 Mark F. Adams <<a href="mailto:mfadams@lbl.gov">mfadams@lbl.gov</a>><br>
> > > > You can try running a problem that finishes with<br>
> > > ><br>
> > > > -malloc<br>
> > > > -malloc_debug<br>
> > > > -malloc_dump<br>
> > > ><br>
> > > > to see if there are any memory issues. It looks like you are<br>
> > completely destroying the solver object each time and starting over.<br>
> > Valgrind would be very useful here because it looks like memory is getting<br>
> > corrupted if the same code runs many times.<br>
> > > ><br>
> > > ><br>
> > > > On Oct 19, 2013, at 3:39 PM, Einar Sørheim <<a href="mailto:einar.sorheim@gmail.com">einar.sorheim@gmail.com</a>><br>
> > wrote:<br>
> > > ><br>
> > > >> yes, I call the given code many times, what i have seen so far is<br>
> > that changing e.g. the order of destruction statements changes the time of<br>
> > failure,I will try to debug to get more details<br>
> > > >><br>
> > > >><br>
> > > >> 2013/10/18 Mark F. Adams <<a href="mailto:mfadams@lbl.gov">mfadams@lbl.gov</a>><br>
> > > >> I would guess there is some memory corruption or leak someplace. Try<br>
> > running with -on_error_attach_debugger, and try to get some more<br>
> > information. Do you just repeat a call to the code below many times? Is<br>
> > this deterministic, does it fail on the same solve each time?<br>
> > > >><br>
> > > >><br>
> > > >> On Oct 18, 2013, at 9:30 AM, Einar Sørheim <<a href="mailto:einar.sorheim@gmail.com">einar.sorheim@gmail.com</a>><br>
> > wrote:<br>
> > > >><br>
> > > >>> Uncommenting KSPSetFromOptions solved my intial problem, however I<br>
> > got a new problem that occurs after many timesteps, in our code we use<br>
> > petsc gamg as a solver for several problems(e.g. thermal and mechanical) at<br>
> > each time step, I have a suspicion it might have something to do with<br>
> > cleanup after solving, the rror appears after many calls to the petsc gamg<br>
> > solver interface, the error message is the following:<br>
> > > >>><br>
> > > >>> [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > > >>> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation<br>
> > Violation, probably memory access out of range<br>
> > > >>> [0]PETSC ERROR: Try option -start_in_debugger or<br>
> > -on_error_attach_debugger<br>
> > > >>> [0]PETSC ERROR: or see<br>
</div></div>> > <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSCERROR" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSCERROR</a>: or try<br>
<div class="HOEnZb"><div class="h5">> > <a href="http://valgrind.org" target="_blank">http://valgrind.org</a> on GNU/linux and Apple Mac OS X to find memory<br>
> > corruption errors<br>
> > > >>> [0]PETSC ERROR: likely location of problem given in stack below<br>
> > > >>> [0]PETSC ERROR: --------------------- Stack Frames<br>
> > ------------------------------------<br>
> > > >>> [0]PETSC ERROR: Note: The EXACT line numbers in the stack are not<br>
> > available,<br>
> > > >>> [0]PETSC ERROR: INSTEAD the line number of the start of the<br>
> > function<br>
> > > >>> [0]PETSC ERROR: is given.<br>
> > > >>> [0]PETSC ERROR: [0] MatSetNearNullSpace line 7580<br>
> > src/mat/interface/D:\PROGRA~1\PETSC-~1.3\src\mat\INTERF~1\matrix.c<br>
> > > >>> [0]PETSC ERROR: --------------------- Error Message<br>
> > ------------------------------------<br>
> > > >>> [0]PETSC ERROR: Signal received!<br>
> > > >>> [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > > >>> [0]PETSC ERROR: Petsc Release Version 3.4.3, Oct, 15, 2013<br>
> > > >>> [0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>
> > > >>> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>
> > > >>> [0]PETSC ERROR: See docs/index.html for manual pages.<br>
> > > >>> [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > > >>> [0]PETSC ERROR: D:\Programming\gits\Alsim\SharedModules -<br>
> > Copy\Source\Release\Alsim_nypetsc.exe on a arch-mswin-c-debug named CMP3 by<br>
> > admeinar Fri Oct 18 12:26:04 2013<br>
> > > >>> [0]PETSC ERROR: Libraries linked from<br>
> > /cygdrive/d/Programming/petsc-3.4.3/arch-mswin-c-debug/lib<br>
> > > >>> [0]PETSC ERROR: Configure run at Thu Oct 17 08:23:18 2013<br>
> > > >>> [0]PETSC ERROR: Configure options --with-debugging=1<br>
> > --with-openmp=yes --with-pthread=no --with-cc="win32fe icl"<br>
> > --with-fc="win32fe ifort" --with-cxx="win32fe icl"<br>
> > --with-blas-lapack-dir="C:/Program Files (x86)/Intel/Composer XE<br>
> > 2013/mkl/lib/intel64" --download-superlu --with-sowing=0 --with-c2html=0<br>
> > --with-mpi-dir="C:/Program Files (x86)/Intel/MPI/<a href="http://4.1.0.028/em64t" target="_blank">4.1.0.028/em64t</a>"<br>
> > --useThreads=0 --useThreads=0<br>
> > > >>><br>
> > > >>> The source code for the petsc interface:<br>
> > > >>><br>
> > > >>> subroutine<br>
> > petscsolver(i_nodes,r_elements_,i_rowstart_,i_columnindex_,r_x_,r_b_,i_Nnz,i_Precond,<br>
> > i_Solver<br>
> > > >>> $ ,coords_,i_nodedof)<br>
> > > >>> #include <finclude/petscsys.h><br>
> > > >>> #include <finclude/petscvec.h><br>
> > > >>> #include <finclude/petscmat.h><br>
> > > >>> #include <finclude/petscpc.h><br>
> > > >>> #include <finclude/petscksp.h><br>
> > > >>> #include <finclude/petscviewer.h><br>
> > > >>><br>
> > > >>> integer i_nodes, i_Nnz, i_Precond, i_solver, i_nodedof<br>
> > > >>> integer i_rowstart_(*),i_columnindex_(*)<br>
> > > >>> double precision r_elements_(*),r_x_(*), r_b_(*),coords_(*)<br>
> > > >>> integer, allocatable :: i_indices_(:)<br>
> > > >>> ! Variables:<br>
> > > >>> !<br>
> > > >>> ! A - matrix that defines linear system<br>
> > > >>> ! ksp - KSP context<br>
> > > >>> ! ksp - KSP context<br>
> > > >>> ! x, b, u - approx solution, RHS, exact solution vectors<br>
> > > >>> !<br>
> > > >>> ! implicit none<br>
> > > >>> Vec x,u,b,vec_coords<br>
> > > >>> PC pc<br>
> > > >>> Mat A,F,Matnull<br>
> > > >>> KSP ksp<br>
> > > >>> PetscInt i,j,II,JJ,m,n,i1,imaxiter<br>
> > > >>> PetscInt Istart,Iend<br>
> > > >>> PetscInt nsteps,one, neg_one<br>
> > > >>> PetscErrorCode ierr<br>
> > > >>> PetscBool flg<br>
> > > >>> PetscScalar v,x_array_(1), norm, tol, t0,t1<br>
> > > >>> PetscOffset i_x<br>
> > > >>> PetscViewer viewer<br>
> > > >>><br>
> > > >>><br>
> > > >>> call PetscTime(t0,ierr)<br>
> > > >>> write(*,*) "Enter petsc"<br>
> > > >>><br>
> > > >>> n = 1<br>
> > > >>> nsteps = 1<br>
> > > >>> one = 1<br>
> > > >>> neg_one = 1<br>
> > > >>> i1 = 1<br>
> > > >>> imaxiter = 1200<br>
> > > >>><br>
> > > >>><br>
> > > >>> call MatCreate(PETSC_COMM_WORLD,A,ierr)<br>
> > > >>> call<br>
> > MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,i_nodes,i_nodes,ierr)<br>
> > > >>> call MatSetType(A, MATSEQAIJ,ierr)<br>
> > > >>> call MatSetUp(A,ierr)<br>
> > > >>> call MatSeqAIJSetPreallocation(A,100,PETSC_NULL_INTEGER,ierr)<br>
> > > >>> ! The matrix is partitioned by contiguous chunks of rows across the<br>
> > > >>> ! processors. Determine which rows of the matrix are locally owned.<br>
> > > >>><br>
> > > >>> call MatGetOwnershipRange(A,Istart,Iend,ierr)<br>
> > > >>><br>
> > > >>><br>
> > > >>> write(*,*) "Start of matrix fill", i_nodes, Istart, Iend<br>
> > > >>> do i=1,i_nodes<br>
> > > >>> ! write(*,*) "Start of loop",i !!, i_indices_(i+1)<br>
> > > >>> ! i_indices_(i+1) = i<br>
> > > >>> ! write(*,*) "Start of loop2"<br>
> > > >>> do j=i_rowstart_(i),i_rowstart_(i+1)-1<br>
> > > >>> ! write(*,*) "Start of loop 3"<br>
> > > >>> v = r_elements_(j)<br>
> > > >>> ! write(*,*) i1,i,i1,j-1,v<br>
> > > >>> call<br>
> > MatSetValues(A,i1,i-1,i1,i_columnindex_(j)-1,v,INSERT_VALUES,ierr)<br>
> > > >>> if (ierr.gt.0) stop<br>
> > > >>> end do<br>
> > > >>> end do<br>
> > > >>> ! write(*,*) "End of matrix fill"<br>
> > > >>><br>
> > > >>> ! Assemble matrix, using the 2-step process:<br>
> > > >>> ! MatAssemblyBegin(), MatAssemblyEnd()<br>
> > > >>> ! Computations can be done while messages are in transition<br>
> > > >>> ! by placing code between these two statements.<br>
> > > >>><br>
> > > >>> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)<br>
> > > >>> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)<br>
> > > >>><br>
> > > >>> ! Create parallel vectors.<br>
> > > >>> ! - When using VecCreate(), the parallel partitioning of the vector<br>
> > > >>> ! is determined by PETSc at runtime.<br>
> > > >>> ! - Note: We form 1 vector from scratch and then duplicate as<br>
> > needed.<br>
> > > >>> ! write(*,*) "Start of vector fill"<br>
> > > >>> call VecCreate(PETSC_COMM_WORLD,u,ierr)<br>
> > > >>> call VecSetType(u, VECSEQ,ierr)<br>
> > > >>> call VecSetSizes(u,PETSC_DECIDE,i_nodes,ierr)<br>
> > > >>> ! call VecSetFromOptions(u,ierr)<br>
> > > >>> call VecDuplicate(u,b,ierr)<br>
> > > >>> call VecDuplicate(b,x,ierr)<br>
> > > >>> do i=1,i_nodes<br>
> > > >>> call VecSetValues(x,i1,i-1,r_x_(i),INSERT_VALUES,ierr)<br>
> > > >>> call VecSetValues(b,i1,i-1,r_b_(i),INSERT_VALUES,ierr)<br>
> > > >>> enddo<br>
> > > >>> ! Assemble vector<br>
> > > >>><br>
> > > >>> call VecAssemblyBegin(x,ierr)<br>
> > > >>> call VecAssemblyEnd(x,ierr)<br>
> > > >>> call VecAssemblyBegin(b,ierr)<br>
> > > >>> call VecAssemblyEnd(b,ierr)<br>
> > > >>> call PetscTime(t1,ierr)<br>
> > > >>> ! write(*,*) "End of vec fill time spent :", t1-t0<br>
> > > >>> ! Create linear solver context<br>
> > > >>><br>
> > > >>> call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)<br>
> > > >>><br>
> > > >>><br>
> > > >>> ! Set operators. Here the matrix that defines the linear system<br>
> > > >>> ! also serves as the preconditioning matrix.<br>
> > > >>><br>
> > > >>> call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr)<br>
> > > >>><br>
> > > >>><br>
> > > >>> call KSPGetPC(ksp,pc,ierr)<br>
> > > >>> tol = 1.e-20<br>
> > > >>> call<br>
> > KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_DOUBLE_PRECISION,imaxiter,ierr)<br>
> > > >>> select case (i_Precond)<br>
> > > >>> case ( 1 )<br>
> > > >>> call PCSetType(pc,PCJACOBI,ierr)<br>
> > > >>> case ( 2 )<br>
> > > >>> call PCSetType(pc,PCILU,ierr)<br>
> > > >>> call PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU,ierr);<br>
> > > >>> call PCFactorSetUpMatSolverPackage(pc,ierr);<br>
> > > >>> call PCFactorGetMatrix(pc,F);<br>
> > > >>> call MatSuperluSetILUDropTol(F,1.e-4);<br>
> > > >>> case ( 3 )<br>
> > > >>> call PCSetType(pc,PCGAMG,ierr)<br>
> > > >>> ! call PCGAMGInitializePackage()<br>
> > > >>> ! call PCCreate_GAMG(pc,ierr)<br>
> > > >>> case DEFAULT<br>
> > > >>> call PCSetType(pc,PCJACOBI,ierr)<br>
> > > >>> end select<br>
> > > >>><br>
> > > >>> ! call PCSetType(pc,PCJACOBI,ierr)<br>
> > > >>> select case (i_solver)<br>
> > > >>> case ( 0 )<br>
> > > >>> call KSPSetType(ksp,KSPCG,ierr)<br>
> > > >>> case DEFAULT<br>
> > > >>> call KSPSetType(ksp,KSPCG,ierr)<br>
> > > >>> end select<br>
> > > >>> call KSPCGSetType(ksp,KSP_CG_SYMMETRIC,ierr)<br>
> > > >>> if (i_nodedof==3) then<br>
> > > >>> write(*,*) "Set coord, number of nodes is:",<br>
> > i_nodes/i_nodedof<br>
> > > >>> call<br>
> > VecCreateSeqWithArray(MPI_COMM_SELF,3,i_nodes,coords_,vec_coords,ierr)<br>
> > > >>> call MatNullSpaceCreateRigidBody(vec_coords,Matnull,ierr)<br>
> > > >>> call MatSetNearNullSpace(A,Matnull,ierr)<br>
> > > >>> call MatNullSpaceDestroy(Matnull,ierr)<br>
> > > >>><br>
> > > >>> write(*,*) "Finish setting null space ierr =", ierr<br>
> > > >>> ! call PCSetCoordinates( pc, 3,i_nodes/i_nodedof, coords_,<br>
> > ierr )<br>
> > > >>><br>
> > > >>> end if<br>
> > > >>> ! call KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL,<br>
> > PETSC_NULL,ierr)<br>
> > > >>> ! Set runtime options, e.g.,<br>
> > > >>> ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol><br>
> > > >>> ! These options will override those specified above as long as<br>
> > > >>> ! KSPSetFromOptions() is called _after_ any other customization<br>
> > > >>> ! routines.<br>
> > > >>><br>
> > > >>> call KSPSetFromOptions(ksp,ierr)<br>
> > > >>> call KSPSetInitialGuessNonzero(ksp,PETSC_TRUE,ierr)<br>
> > > >>> call KSPSetUp(ksp,ierr)<br>
> > > >>> call PetscTime(t0,ierr)<br>
> > > >>> write(*,*) "Time setup", t0-t1<br>
> > > >>> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> > - -<br>
> > > >>> ! Solve the linear system<br>
> > > >>> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
> > - -<br>
> > > >>><br>
> > > >>> call KSPSolve(ksp,b,x,ierr)<br>
> > > >>> call PetscTime(t1,ierr)<br>
> > > >>> write(*,*) "end solve, time used:",t1-t0<br>
> > > >>> call KSPGetSolution(ksp,x,ierr)<br>
> > > >>> call VecGetArray(x,x_array_,i_x,ierr)<br>
> > > >>> do i=1,i_nodes<br>
> > > >>> r_x_(i) = x_array_(i_x + i)<br>
> > > >>> enddo<br>
> > > >>> call VecRestoreArray(x,x_array_,i_x,ierr)<br>
> > > >>><br>
> > > >>><br>
> > > >>> ! Free work space. All PETSc objects should be destroyed when they<br>
> > > >>> ! are no longer needed.<br>
> > > >>> ! call PCDestroy(pc,ierr)<br>
> > > >>><br>
> > > >>> call VecDestroy(u,ierr)<br>
> > > >>> call VecDestroy(x,ierr)<br>
> > > >>> call VecDestroy(b,ierr)<br>
> > > >>> call MatDestroy(A,ierr)<br>
> > > >>> if (i_nodedof==3) then<br>
> > > >>> call VecDestroy(vec_coords,ierr)<br>
> > > >>> call MatDestroy(Matnull,ierr)<br>
> > > >>> endif<br>
> > > >>> call KSPDestroy(ksp,ierr)<br>
> > > >>> call PetscTime(t0,ierr)<br>
> > > >>> write(*,*) "time celan up :", t0-t1<br>
> > > >>><br>
> > > >>> end<br>
> > > >>><br>
> > > >>> --<br>
> > > >>> Einar Sørheim<br>
> > > >>> <Mail Attachment><br>
> > > >><br>
> > > >><br>
> > > >><br>
> > > >><br>
> > > >> --<br>
> > > >> Einar Sørheim<br>
> > > ><br>
> > > ><br>
> > > ><br>
> > > ><br>
> > > > --<br>
> > > > Einar Sørheim<br>
> > ><br>
> > ><br>
> > ><br>
> > ><br>
> > > --<br>
> > > Einar Sørheim<br>
> > > <log_summary.txt><br>
> ><br>
> ><br>
><br>
><br>
><br>
</div></div></blockquote></div><br></div>