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