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