! ! "$Id: ex6f.F,v 1.36 2001/08/07 03:04:00 balay Exp $"; ! ! Description: This example demonstrates repeated linear solves as ! well as the use of different preconditioner and linear system ! matrices. This example also illustrates how to save PETSc objects ! in common blocks. ! !/*T ! Concepts: SLES^repeatedly solving linear systems; ! Concepts: SLES^different matrices for linear system and preconditioner; ! Processors: n !T*/ ! ! The following include statements are required for SLES Fortran programs: ! petsc.h - base PETSc routines ! petscvec.h - vectors ! petscmat.h - matrices ! petscpc.h - preconditioners ! petscksp.h - Krylov subspace methods ! petscsles.h - SLES interface ! Other include statements may be needed if using additional PETSc ! routines in a Fortran program, e.g., ! petscviewer.h - viewers ! petscis.h - index sets ! program main #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscmat.h" #include "include/finclude/petscsles.h" #include "include/finclude/petscpc.h" #include "include/finclude/petscksp.h" ! Variables: ! ! A - matrix that defines linear system ! sles - SLES context ! ksp - KSP context ! x, b, u - approx solution, RHS, exact solution vectors ! Vec x,u,b Mat A SLES sles integer i,j,II,JJ,ierr,m,n integer Istart,Iend,flg,nsteps PetscScalar v call PetscInitialize(PETSC_NULL_CHARACTER,ierr) m = 3 n = 3 nsteps = 2 call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr) call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr) call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-nsteps',nsteps, & & flg,ierr) ! Create parallel matrix, specifying only its global dimensions. ! When using MatCreate(), the matrix format can be specified at ! runtime. Also, the parallel partitioning of the matrix is ! determined by PETSc at runtime. call MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m*n, & & m*n,A,ierr) call MatSetFromOptions(A,ierr) ! The matrix is partitioned by contiguous chunks of rows across the ! processors. Determine which rows of the matrix are locally owned. call MatGetOwnershipRange(A,Istart,Iend,ierr) ! Set matrix elements. ! - Each processor needs to insert only elements that it owns ! locally (but any non-local elements will be sent to the ! appropriate processor during matrix assembly). ! - Always specify global rows and columns of matrix entries. do 10, II=Istart,Iend-1 v = -1.0 i = II/n j = II - i*n if (i.gt.0) then JJ = II - n call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) endif if (i.lt.m-1) then JJ = II + n call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) endif if (j.gt.0) then JJ = II - 1 call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) endif if (j.lt.n-1) then JJ = II + 1 call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) endif v = 4.0 call MatSetValues(A,1,II,1,II,v,ADD_VALUES,ierr) 10 continue ! Assemble matrix, using the 2-step process: ! MatAssemblyBegin(), MatAssemblyEnd() ! Computations can be done while messages are in transition ! by placing code between these two statements. call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) ! Create parallel vectors. ! - When using VecCreate(), the parallel partitioning of the vector ! is determined by PETSc at runtime. ! - Note: We form 1 vector from scratch and then duplicate as needed. call VecCreate(PETSC_COMM_WORLD,u,ierr) call VecSetSizes(u,PETSC_DECIDE,m*n,ierr) call VecSetFromOptions(u,ierr) call VecDuplicate(u,b,ierr) call VecDuplicate(b,x,ierr) ! Create linear solver context call SLESCreate(PETSC_COMM_WORLD,sles,ierr) ! Set runtime options (e.g., -ksp_type -pc_type ) call SLESSetFromOptions(sles,ierr) ! Solve several linear systems in succession do 100 i=1,nsteps call solve1(sles,A,x,b,u,i,nsteps,ierr) 100 continue ! Free work space. All PETSc objects should be destroyed when they ! are no longer needed. call VecDestroy(u,ierr) call VecDestroy(x,ierr) call VecDestroy(b,ierr) call MatDestroy(A,ierr) call SLESDestroy(sles,ierr) call PetscFinalize(ierr) end ! ----------------------------------------------------------------------- ! subroutine solve1(sles,A,x,b,u,count,nsteps,ierr) #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscmat.h" #include "include/finclude/petscsles.h" #include "include/finclude/petscpc.h" #include "include/finclude/petscksp.h" ! ! solve1 - This routine is used for repeated linear system solves. ! We update the linear system matrix each time, but retain the same ! preconditioning matrix for all linear solves. ! ! A - linear system matrix ! A2 - preconditioning matrix ! PetscScalar v,val integer II,ierr,Istart,Iend,count,nsteps,its Mat A SLES sles KSP ksp Vec x,b,u ! Use common block to retain matrix between successive subroutine calls Mat A2 integer rank,pflag common /my_data/ A2,pflag,rank ! First time thorough: Create new matrix to define the linear system if (count .eq. 1) then call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) pflag = 0 call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-mat_view', & & pflag,ierr) if (pflag .ne. 0) then if (rank .eq. 0) write(6,100) endif call MatConvert(A,MATSAME,A2,ierr) ! All other times: Set previous solution as initial guess for next solve. else call SLESGetKSP(sles,ksp,ierr) call KSPSetInitialGuessNonzero(ksp,PETSC_TRUE,ierr) endif ! Alter the matrix A a bit call MatGetOwnershipRange(A,Istart,Iend,ierr) do 20, II=Istart,Iend-1 v = 2.0 call MatSetValues(A,1,II,1,II,v,ADD_VALUES,ierr) 20 continue call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) if (pflag .ne. 0) then if (rank .eq. 0) write(6,110) endif call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) ! Set the exact solution; compute the right-hand-side vector val = 1.0*count call VecSet(val,u,ierr) call MatMult(A,u,b,ierr) ! Set operators, keeping the identical preconditioner matrix for ! all linear solves. This approach is often effective when the ! linear systems do not change very much between successive steps. call SLESSetOperators(sles,A,A2,SAME_PRECONDITIONER,ierr) ! Solve linear system call SLESSolve(sles,b,x,its,ierr) ! Destroy the preconditioner matrix on the last time through if (count .eq. nsteps) call MatDestroy(A2,ierr) 100 format('previous matrix: preconditioning') 110 format('next matrix: defines linear system') end