! ! 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: KSP^repeatedly solving linear systems; ! Concepts: KSP^different matrices for linear system and preconditioner; ! Processors: n !T*/ ! ! The following include statements are required for KSP Fortran programs: ! petscsys.h - base PETSc routines ! petscvec.h - vectors ! petscmat.h - matrices ! petscpc.h - preconditioners ! petscksp.h - Krylov subspace methods ! 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 #include #include #include ! Variables: ! ! A - matrix that defines linear system ! ksp - KSP context ! ksp - KSP context ! x, b, u - approx solution, RHS, exact solution vectors ! Vec x,u,b Mat A KSP ksp PetscInt i,j,II,JJ,m,n PetscInt Istart,Iend PetscInt nsteps,one PetscErrorCode ierr PetscBool flg PetscScalar v call PetscInitialize(PETSC_NULL_CHARACTER,ierr) m = 3 n = 3 nsteps = 2 one = 1 call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & & '-m',m,flg,ierr) call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & & '-n',n,flg,ierr) call PetscOptionsGetInt(PETSC_NULL_OBJECT,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,A,ierr) call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr) call MatSetFromOptions(A,ierr) call MatSetUp(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,one,II,one,JJ,v,ADD_VALUES,ierr) endif if (i.lt.m-1) then JJ = II + n call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) endif if (j.gt.0) then JJ = II - 1 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) endif if (j.lt.n-1) then JJ = II + 1 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) endif v = 4.0 call MatSetValues(A,one,II,one,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 KSPCreate(PETSC_COMM_WORLD,ksp,ierr) ! Set runtime options (e.g., -ksp_type -pc_type ) call KSPSetFromOptions(ksp,ierr) ! Solve several linear systems in succession do 100 i=1,nsteps call solve1(ksp,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 KSPDestroy(ksp,ierr) call PetscFinalize(ierr) end ! ----------------------------------------------------------------------- ! subroutine solve1(ksp,A,x,b,u,count,nsteps,ierr) #include #include #include #include #include ! ! 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 PetscInt II,Istart,Iend PetscInt count,nsteps,one PetscErrorCode ierr Mat A KSP ksp Vec x,b,u ! Use common block to retain matrix between successive subroutine calls Mat A2 PetscMPIInt rank PetscBool pflag common /my_data/ A2,pflag,rank one = 1 ! 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 = .false. call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & & '-mat_view',pflag,ierr) if (pflag) then if (rank .eq. 0) write(6,100) call flush(6) endif call MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,A2,ierr) ! All other times: Set previous solution as initial guess for next solve. else 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,one,II,one,II,v,ADD_VALUES,ierr) 20 continue call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) if (pflag) then if (rank .eq. 0) write(6,110) call flush(6) endif call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) ! Set the exact solution; compute the right-hand-side vector val = 1.0*real(count) call VecSet(u,val,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 KSPSetReusePreconditioner(ksp,PETSC_TRUE,ierr) call KSPSetOperators(ksp,A,A2,ierr) ! Solve linear system call KSPSolve(ksp,b,x,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