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