1*c4762a1bSJed Brown! 2*c4762a1bSJed Brown! Description: Tests PCFieldSplitGetIS and PCFieldSplitSetIS from Fortan. 3*c4762a1bSJed Brown! 4*c4762a1bSJed Brown program main 5*c4762a1bSJed Brown#include <petsc/finclude/petscksp.h> 6*c4762a1bSJed Brown use petscksp 7*c4762a1bSJed Brown implicit none 8*c4762a1bSJed Brown 9*c4762a1bSJed Brown Vec x,b,u 10*c4762a1bSJed Brown Mat A 11*c4762a1bSJed Brown KSP ksp 12*c4762a1bSJed Brown PC pc 13*c4762a1bSJed Brown PetscReal norm; 14*c4762a1bSJed Brown PetscErrorCode ierr 15*c4762a1bSJed Brown PetscInt i,n,col(3),its,i1,i2,i3 16*c4762a1bSJed Brown PetscInt ione,izero 17*c4762a1bSJed Brown PetscBool flg 18*c4762a1bSJed Brown PetscMPIInt size 19*c4762a1bSJed Brown PetscScalar none,one,value(3) 20*c4762a1bSJed Brown IS isin,isout 21*c4762a1bSJed Brown 22*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 23*c4762a1bSJed Brown! Beginning of program 24*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 25*c4762a1bSJed Brown 26*c4762a1bSJed Brown call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 27*c4762a1bSJed Brown if (ierr .ne. 0) then 28*c4762a1bSJed Brown print*,'Unable to initialize PETSc' 29*c4762a1bSJed Brown stop 30*c4762a1bSJed Brown endif 31*c4762a1bSJed Brown call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) 32*c4762a1bSJed Brown if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif 33*c4762a1bSJed Brown none = -1.0 34*c4762a1bSJed Brown one = 1.0 35*c4762a1bSJed Brown n = 10 36*c4762a1bSJed Brown i1 = 1 37*c4762a1bSJed Brown i2 = 2 38*c4762a1bSJed Brown i3 = 3 39*c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr) 40*c4762a1bSJed Brown 41*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 42*c4762a1bSJed Brown! Compute the matrix and right-hand-side vector that define 43*c4762a1bSJed Brown! the linear system, Ax = b. 44*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 45*c4762a1bSJed Brown 46*c4762a1bSJed Brown! Create matrix. When using MatCreate(), the matrix format can 47*c4762a1bSJed Brown! be specified at runtime. 48*c4762a1bSJed Brown 49*c4762a1bSJed Brown call MatCreate(PETSC_COMM_WORLD,A,ierr) 50*c4762a1bSJed Brown call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr) 51*c4762a1bSJed Brown call MatSetFromOptions(A,ierr) 52*c4762a1bSJed Brown call MatSetUp(A,ierr) 53*c4762a1bSJed Brown 54*c4762a1bSJed Brown! Assemble matrix. 55*c4762a1bSJed Brown! - Note that MatSetValues() uses 0-based row and column numbers 56*c4762a1bSJed Brown! in Fortran as well as in C (as set here in the array "col"). 57*c4762a1bSJed Brown 58*c4762a1bSJed Brown value(1) = -1.0 59*c4762a1bSJed Brown value(2) = 2.0 60*c4762a1bSJed Brown value(3) = -1.0 61*c4762a1bSJed Brown do 50 i=1,n-2 62*c4762a1bSJed Brown col(1) = i-1 63*c4762a1bSJed Brown col(2) = i 64*c4762a1bSJed Brown col(3) = i+1 65*c4762a1bSJed Brown call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr) 66*c4762a1bSJed Brown 50 continue 67*c4762a1bSJed Brown i = n - 1 68*c4762a1bSJed Brown col(1) = n - 2 69*c4762a1bSJed Brown col(2) = n - 1 70*c4762a1bSJed Brown call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr) 71*c4762a1bSJed Brown i = 0 72*c4762a1bSJed Brown col(1) = 0 73*c4762a1bSJed Brown col(2) = 1 74*c4762a1bSJed Brown value(1) = 2.0 75*c4762a1bSJed Brown value(2) = -1.0 76*c4762a1bSJed Brown call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr) 77*c4762a1bSJed Brown call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) 78*c4762a1bSJed Brown call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) 79*c4762a1bSJed Brown 80*c4762a1bSJed Brown! Create vectors. Note that we form 1 vector from scratch and 81*c4762a1bSJed Brown! then duplicate as needed. 82*c4762a1bSJed Brown 83*c4762a1bSJed Brown call VecCreate(PETSC_COMM_WORLD,x,ierr) 84*c4762a1bSJed Brown call VecSetSizes(x,PETSC_DECIDE,n,ierr) 85*c4762a1bSJed Brown call VecSetFromOptions(x,ierr) 86*c4762a1bSJed Brown call VecDuplicate(x,b,ierr) 87*c4762a1bSJed Brown call VecDuplicate(x,u,ierr) 88*c4762a1bSJed Brown 89*c4762a1bSJed Brown! Set exact solution; then compute right-hand-side vector. 90*c4762a1bSJed Brown 91*c4762a1bSJed Brown call VecSet(u,one,ierr) 92*c4762a1bSJed Brown call MatMult(A,u,b,ierr) 93*c4762a1bSJed Brown 94*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95*c4762a1bSJed Brown! Create the linear solver and set various options 96*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97*c4762a1bSJed Brown 98*c4762a1bSJed Brown! Create linear solver context 99*c4762a1bSJed Brown 100*c4762a1bSJed Brown call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) 101*c4762a1bSJed Brown 102*c4762a1bSJed Brown! Set operators. Here the matrix that defines the linear system 103*c4762a1bSJed Brown! also serves as the preconditioning matrix. 104*c4762a1bSJed Brown 105*c4762a1bSJed Brown call KSPSetOperators(ksp,A,A,ierr) 106*c4762a1bSJed Brown 107*c4762a1bSJed Brown! Set linear solver defaults for this problem (optional). 108*c4762a1bSJed Brown! - By extracting the KSP and PC contexts from the KSP context, 109*c4762a1bSJed Brown! we can then directly directly call any KSP and PC routines 110*c4762a1bSJed Brown! to set various options. 111*c4762a1bSJed Brown! - The following four statements are optional; all of these 112*c4762a1bSJed Brown! parameters could alternatively be specified at runtime via 113*c4762a1bSJed Brown! KSPSetFromOptions(); 114*c4762a1bSJed Brown 115*c4762a1bSJed Brown call KSPGetPC(ksp,pc,ierr) 116*c4762a1bSJed Brown call PCSetType(pc,PCFIELDSPLIT,ierr) 117*c4762a1bSJed Brown izero = 0 118*c4762a1bSJed Brown ione = 1 119*c4762a1bSJed Brown call ISCreateStride(PETSC_COMM_SELF,n,izero,ione,isin,ierr) 120*c4762a1bSJed Brown call PCFieldSplitSetIS(pc,"splitname",isin,ierr) 121*c4762a1bSJed Brown call PCFieldSplitGetIS(pc,"splitname",isout,ierr) 122*c4762a1bSJed Brown if (isin .ne. isout) then ; SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCFieldSplitGetIS() failed"); endif 123*c4762a1bSJed Brown 124*c4762a1bSJed Brown! Set runtime options, e.g., 125*c4762a1bSJed Brown! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> 126*c4762a1bSJed Brown! These options will override those specified above as long as 127*c4762a1bSJed Brown! KSPSetFromOptions() is called _after_ any other customization 128*c4762a1bSJed Brown! routines. 129*c4762a1bSJed Brown 130*c4762a1bSJed Brown call KSPSetFromOptions(ksp,ierr) 131*c4762a1bSJed Brown 132*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 133*c4762a1bSJed Brown! Solve the linear system 134*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135*c4762a1bSJed Brown 136*c4762a1bSJed Brown call KSPSolve(ksp,b,x,ierr) 137*c4762a1bSJed Brown call PetscLogStagePop(ierr) 138*c4762a1bSJed Brown 139*c4762a1bSJed Brown! View solver info; we could instead use the option -ksp_view 140*c4762a1bSJed Brown 141*c4762a1bSJed Brown call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr) 142*c4762a1bSJed Brown 143*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 144*c4762a1bSJed Brown! Check solution and clean up 145*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146*c4762a1bSJed Brown 147*c4762a1bSJed Brown! Check the error 148*c4762a1bSJed Brown 149*c4762a1bSJed Brown call VecAXPY(x,none,u,ierr) 150*c4762a1bSJed Brown call VecNorm(x,NORM_2,norm,ierr) 151*c4762a1bSJed Brown call KSPGetIterationNumber(ksp,its,ierr) 152*c4762a1bSJed Brown if (norm .gt. 1.e-12) then 153*c4762a1bSJed Brown write(6,100) norm,its 154*c4762a1bSJed Brown else 155*c4762a1bSJed Brown write(6,200) its 156*c4762a1bSJed Brown endif 157*c4762a1bSJed Brown 100 format('Norm of error ',e11.4,', Iterations = ',i5) 158*c4762a1bSJed Brown 200 format('Norm of error < 1.e-12, Iterations = ',i5) 159*c4762a1bSJed Brown 160*c4762a1bSJed Brown! Free work space. All PETSc objects should be destroyed when they 161*c4762a1bSJed Brown! are no longer needed. 162*c4762a1bSJed Brown 163*c4762a1bSJed Brown call ISDestroy(isin,ierr) 164*c4762a1bSJed Brown call VecDestroy(x,ierr) 165*c4762a1bSJed Brown call VecDestroy(u,ierr) 166*c4762a1bSJed Brown call VecDestroy(b,ierr) 167*c4762a1bSJed Brown call MatDestroy(A,ierr) 168*c4762a1bSJed Brown call KSPDestroy(ksp,ierr) 169*c4762a1bSJed Brown call PetscFinalize(ierr) 170*c4762a1bSJed Brown 171*c4762a1bSJed Brown end 172*c4762a1bSJed Brown 173*c4762a1bSJed Brown!/*TEST 174*c4762a1bSJed Brown! 175*c4762a1bSJed Brown! test: 176*c4762a1bSJed Brown! args: -ksp_monitor 177*c4762a1bSJed Brown! 178*c4762a1bSJed Brown!TEST*/ 179