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