program main implicit none #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscmat.h" #include "include/finclude/petscpc.h" #include "include/finclude/petscksp.h" #include "include/finclude/petscsys.h" ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Variable declarations ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! Variables: ! ksp - linear solver context ! ksp - Krylov subspace method context ! pc - preconditioner context ! x, b, u - approx solution, right-hand-side, exact solution vectors ! A - matrix that defines linear system ! its - iterations for convergence ! norm - norm of error in solution ! rctx - random number generator context ! double precision norm PetscInt :: i,j,ii,jj,m,n,its PetscInt :: istart,iend,ione PetscErrorCode :: ierr PetscMPIInt :: rank, size PetscTruth :: flg PetscScalar :: v,one,neg_one Vec :: es,hs,hb,eb Mat :: z1,z2 KSP :: eksp,hksp PetscRandom :: rctx external kspmon, kspconverged ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Beginning of program ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Initsection ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call petscinitialize(petsc_null_character,ierr) m = 3 n = 3 one = 1.0 neg_one = -1.0 ione = 1 call petscoptionsgetint(petsc_null_character,'-m',m,flg,ierr) call petscoptionsgetint(petsc_null_character,'-n',n,flg,ierr) call mpi_comm_rank(petsc_comm_world,rank,ierr) call mpi_comm_size(petsc_comm_world,size,ierr) call matcreate(petsc_comm_world,z1,ierr) call matcreate(petsc_comm_world,z2,ierr) call matsetsizes(z1,petsc_decide,petsc_decide,m*n,m*n,ierr) call matsetsizes(z2,petsc_decide,petsc_decide,m*n,m*n,ierr) call matsetfromoptions(z1,ierr) call matsetfromoptions(z2,ierr) call matgetownershiprange(z1,istart,iend,ierr) call matgetownershiprange(z2,istart,iend,ierr) do ii=istart,iend-1 v = -1.0 i = ii/n j = ii - i*n if (i.gt.0) then jj = ii - n call matsetvalues(z1,ione,ii,ione,jj,v,insert_values,ierr) endif if (i.lt.m-1) then jj = ii + n call matsetvalues(z1,ione,ii,ione,jj,v,insert_values,ierr) endif if (j.gt.0) then jj = ii - 1 call matsetvalues(z1,ione,ii,ione,jj,v,insert_values,ierr) endif if (j.lt.n-1) then jj = ii + 1 call matsetvalues(z1,ione,ii,ione,jj,v,insert_values,ierr) endif v = 4.0 call matsetvalues(z1,ione,ii,ione,ii,v,insert_values,ierr) end do call matassemblybegin(z1,mat_final_assembly,ierr) do ii=istart,iend-1 v = -1.0 i = ii/n j = ii - i*n if (i.gt.0) then jj = ii - n call matsetvalues(z2,ione,ii,ione,jj,v,insert_values,ierr) endif if (i.lt.m-1) then jj = ii + n call matsetvalues(z2,ione,ii,ione,jj,v,insert_values,ierr) endif if (j.gt.0) then jj = ii - 1 call matsetvalues(z2,ione,ii,ione,jj,v,insert_values,ierr) endif if (j.lt.n-1) then jj = ii + 1 call matsetvalues(z2,ione,ii,ione,jj,v,insert_values,ierr) endif v = 4.0 call matsetvalues(z2,ione,ii,ione,ii,v,insert_values,ierr) end do ! your code can go here. call matassemblybegin(z2,mat_final_assembly,ierr) call veccreatempi(petsc_comm_world,petsc_decide,m*n,hb,ierr) call vecsetfromoptions(hb,ierr) call vecduplicate(hb,es,ierr) call vecduplicate(es,eb,ierr) call vecduplicate(eb,hs,ierr) call matassemblyend(z1,mat_final_assembly,ierr) call matassemblyend(z2,mat_final_assembly,ierr) !call PetscOptionsHasName(PETSC_NULL_CHARACTER,"-random_exact_sol",flg,ierr) !if (flg .eq. 1) then ! call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr) ! call PetscRandomSetFromOptions(rctx,ierr) ! call VecSetRandom(u,rctx,ierr) ! call PetscRandomDestroy(rctx,ierr) !else call vecset(hb,one,ierr) !endif !call MatMult(A,u,b,ierr) ! View the exact solution vector if desired !call PetscOptionsHasName(PETSC_NULL_CHARACTER,"-view_exact_sol",flg,ierr) !if (flg .eq. 1) then ! call VecView(u,PETSC_VIEWER_STDOUT_WORLD,ierr) !endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create the linear solver and set various options ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create linear solver context call kspcreate(petsc_comm_world,eksp,ierr) call kspcreate(petsc_comm_world,hksp,ierr) ! Set operators. Here the matrix that defines the linear system ! also serves as the preconditioning matrix. ! call KSPGetPC(ksp,pc,ierr) ! ptype = PCJACOBI ! call PCSetType(pc,ptype,ierr) ! tol = 1.e-7 ! call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, & ! & PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr) ! Set user-defined monitoring routine if desired !call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-my_ksp_monitor',flg,ierr) !if (flg .eq. 1) then ! call KSPMonitorSet(ksp,MyKSPMonitor,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) !endif call kspsetfromoptions(eksp,ierr) call kspsetfromoptions(hksp,ierr) ! Set convergence test routine if desired !call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-my_ksp_convergence',flg,ierr) !if (flg .eq. 1) then ! call KSPSetConvergenceTest(ksp,MyKSPConverged,PETSC_NULL_OBJECT,ierr) !endif ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Solve the linear system ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - do ii=1,10 call kspsetoperators(eksp,z1,z1,different_nonzero_pattern,ierr) call kspsolve(eksp,eb,es,ierr) ! ! Code missing! ! call kspsetoperators(hksp,z2,z2,different_nonzero_pattern,ierr) call kspsolve(hksp,hb,hs,ierr) end do ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Check solution and clean up ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Check the error !call VecAXPY(x,neg_one,u,ierr) !call VecNorm(x,NORM_2,norm,ierr) !call KSPGetIterationNumber(ksp,its,ierr) !if (rank .eq. 0) then ! if (norm .gt. 1.e-12) then ! write(6,100) norm,its ! else ! write(6,110) its ! endif !endif !100 format('Norm of error ',e10.4,' iterations ',i5) !110 format('Norm of error < 1.e-12,iterations ',i5) ! Sleep for a while so you can check the graphs !call matview(z1,PETSC_VIEWER_DRAW_WORLD,ierr) !call petscsleep(10,ierr) call kspdestroy(eksp,ierr) call kspdestroy(hksp,ierr) call vecdestroy(es,ierr) call vecdestroy(hs,ierr) call vecdestroy(eb,ierr) call vecdestroy(hb,ierr) call matdestroy(z1,ierr) call matdestroy(z2,ierr) ! Always call PetscFinalize() before exiting a program. This routine call petscfinalize(ierr) end program main ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Subroutines below here ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine kspmon(ksp,n,rnorm,dummy,ierr) implicit none #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscksp.h" KSP :: ksp Vec :: x PetscErrorCode :: ierr PetscInt :: n,dummy PetscMPIInt :: rank double precision :: rnorm ! Build the solution vector call kspbuildsolution(ksp,petsc_null_object,x,ierr) ! Write the solution vector and residual norm to stdout ! - Note that the parallel viewer PETSC_VIEWER_STDOUT_WORLD ! handles data from multiple processors so that the ! output is not jumbled. call mpi_comm_rank(petsc_comm_world,rank,ierr) if (rank .eq. 0) write(6,100) n call vecview(x,petsc_viewer_stdout_world,ierr) if (rank .eq. 0) write(6,200) n,rnorm 100 format('iteration ',i5,' solution vector:') 200 format('iteration ',i5,' residual norm ',e10.4) ierr = 0 end subroutine kspmon ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine kspconverged(ksp,n,rnorm,flag,dummy,ierr) implicit none #include "include/finclude/petsc.h" #include "include/finclude/petscvec.h" #include "include/finclude/petscksp.h" KSP :: ksp PetscErrorCode :: ierr PetscInt :: n,dummy KSPConvergedReason :: flag double precision :: rnorm if (rnorm .le. .05) then flag = 1 else flag = 0 endif ierr = 0 end subroutine kspconverged ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Below is nothing ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -