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