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 134! View solver info; we could instead use the option -ksp_view 135 136 PetscCallA(KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)) 137 138! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139! Check solution and clean up 140! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 142! Check the error 143 144 PetscCallA(VecAXPY(x,none,u,ierr)) 145 PetscCallA(VecNorm(x,NORM_2,norm,ierr)) 146 PetscCallA(KSPGetIterationNumber(ksp,its,ierr)) 147 if (norm .gt. 1.e-12) then 148 write(6,100) norm,its 149 else 150 write(6,200) its 151 endif 152 100 format('Norm of error ',e11.4,', Iterations = ',i5) 153 200 format('Norm of error < 1.e-12, Iterations = ',i5) 154 155! Free work space. All PETSc objects should be destroyed when they 156! are no longer needed. 157 158 PetscCallA(ISDestroy(isin,ierr)) 159 PetscCallA(VecDestroy(x,ierr)) 160 PetscCallA(VecDestroy(u,ierr)) 161 PetscCallA(VecDestroy(b,ierr)) 162 PetscCallA(MatDestroy(A,ierr)) 163 PetscCallA(KSPDestroy(ksp,ierr)) 164 PetscCallA(PetscFinalize(ierr)) 165 166 end 167 168!/*TEST 169! 170! test: 171! args: -ksp_monitor 172! 173!TEST*/ 174