1 #include <petsctao.h> /*I "petsctao.h" I*/ 2 #include <petsc-private/matimpl.h> 3 #include <../src/tao/matrix/submatfree.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "VecGetSubVec" 7 /*@ 8 VecGetSubVec - Gets a subvector using the IS 9 10 Input Parameters: 11 + vfull - the full matrix 12 . is - the index set for the subvector 13 . reduced_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, TAO_SUBSET_MATRIXFREE) 14 - maskvalue - the value to set the unused vector elements to (for TAO_SUBSET_MASK or TAO_SUBSET_MATRIXFREE) 15 16 Output Parameters: 17 . vreduced - the subvector 18 19 Note: 20 maskvalue should usually be 0.0, unless a pointwise divide will be used. 21 @*/ 22 PetscErrorCode VecGetSubVec(Vec vfull, IS is, PetscInt reduced_type, PetscReal maskvalue, Vec *vreduced) 23 { 24 PetscErrorCode ierr; 25 PetscInt nfull,nreduced,nreduced_local,rlow,rhigh,flow,fhigh; 26 PetscInt i,nlocal; 27 PetscReal *fv,*rv; 28 const PetscInt *s; 29 IS ident; 30 VecType vtype; 31 VecScatter scatter; 32 MPI_Comm comm; 33 34 PetscFunctionBegin; 35 PetscValidHeaderSpecific(vfull,VEC_CLASSID,1); 36 PetscValidHeaderSpecific(is,IS_CLASSID,2); 37 38 ierr = VecGetSize(vfull, &nfull);CHKERRQ(ierr); 39 ierr = ISGetSize(is, &nreduced);CHKERRQ(ierr); 40 41 if (nreduced == nfull) { 42 ierr = VecDestroy(vreduced);CHKERRQ(ierr); 43 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr); 44 ierr = VecCopy(vfull,*vreduced);CHKERRQ(ierr); 45 } else { 46 switch (reduced_type) { 47 case TAO_SUBSET_SUBVEC: 48 ierr = VecGetType(vfull,&vtype);CHKERRQ(ierr); 49 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr); 50 ierr = ISGetLocalSize(is,&nreduced_local);CHKERRQ(ierr); 51 ierr = PetscObjectGetComm((PetscObject)vfull,&comm);CHKERRQ(ierr); 52 if (*vreduced) { 53 ierr = VecDestroy(vreduced);CHKERRQ(ierr); 54 } 55 ierr = VecCreate(comm,vreduced);CHKERRQ(ierr); 56 ierr = VecSetType(*vreduced,vtype);CHKERRQ(ierr); 57 58 ierr = VecSetSizes(*vreduced,nreduced_local,nreduced);CHKERRQ(ierr); 59 ierr = VecGetOwnershipRange(*vreduced,&rlow,&rhigh);CHKERRQ(ierr); 60 ierr = ISCreateStride(comm,nreduced_local,rlow,1,&ident);CHKERRQ(ierr); 61 ierr = VecScatterCreate(vfull,is,*vreduced,ident,&scatter);CHKERRQ(ierr); 62 ierr = VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 63 ierr = VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 64 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 65 ierr = ISDestroy(&ident);CHKERRQ(ierr); 66 break; 67 68 case TAO_SUBSET_MASK: 69 case TAO_SUBSET_MATRIXFREE: 70 /* vr[i] = vf[i] if i in is 71 vr[i] = 0 otherwise */ 72 if (*vreduced == NULL) { 73 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr); 74 } 75 76 ierr = VecSet(*vreduced,maskvalue);CHKERRQ(ierr); 77 ierr = ISGetLocalSize(is,&nlocal);CHKERRQ(ierr); 78 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr); 79 ierr = VecGetArray(vfull,&fv);CHKERRQ(ierr); 80 ierr = VecGetArray(*vreduced,&rv);CHKERRQ(ierr); 81 ierr = ISGetIndices(is,&s);CHKERRQ(ierr); 82 if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_WORLD,1,"IS local size %d > Vec local size %d",nlocal,fhigh-flow); 83 for (i=0;i<nlocal;i++) { 84 rv[s[i]-flow] = fv[s[i]-flow]; 85 } 86 ierr = ISRestoreIndices(is,&s);CHKERRQ(ierr); 87 ierr = VecRestoreArray(vfull,&fv);CHKERRQ(ierr); 88 ierr = VecRestoreArray(*vreduced,&rv);CHKERRQ(ierr); 89 break; 90 } 91 } 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "MatGetSubMat" 97 /*@ 98 MatGetSubMat - Gets a submatrix using the IS 99 100 Input Parameters: 101 + M - the full matrix (n x n) 102 . is - the index set for the submatrix (both row and column index sets need to be the same) 103 . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option 104 - subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, 105 TAO_SUBSET_MATRIXFREE) 106 107 Output Parameters: 108 . Msub - the submatrix 109 @*/ 110 PetscErrorCode MatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub) 111 { 112 PetscErrorCode ierr; 113 IS iscomp; 114 PetscBool flg; 115 116 PetscFunctionBegin; 117 PetscValidHeaderSpecific(M,MAT_CLASSID,1); 118 PetscValidHeaderSpecific(is,IS_CLASSID,2); 119 ierr = MatDestroy(Msub);CHKERRQ(ierr); 120 switch (subset_type) { 121 case TAO_SUBSET_SUBVEC: 122 ierr = MatGetSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);CHKERRQ(ierr); 123 break; 124 125 case TAO_SUBSET_MASK: 126 /* Get Reduced Hessian 127 Msub[i,j] = M[i,j] if i,j in Free_Local or i==j 128 Msub[i,j] = 0 if i!=j and i or j not in Free_Local 129 */ 130 ierr = PetscOptionsBool("-different_submatrix","use separate hessian matrix when computing submatrices","TaoSubsetType",PETSC_FALSE,&flg,NULL); 131 if (flg == PETSC_TRUE) { 132 ierr = MatDuplicate(M, MAT_COPY_VALUES, Msub);CHKERRQ(ierr); 133 } else { 134 /* Act on hessian directly (default) */ 135 ierr = PetscObjectReference((PetscObject)M);CHKERRQ(ierr); 136 *Msub = M; 137 } 138 /* Save the diagonal to temporary vector */ 139 ierr = MatGetDiagonal(*Msub,v1);CHKERRQ(ierr); 140 141 /* Zero out rows and columns */ 142 ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr); 143 144 /* Use v1 instead of 0 here because of PETSc bug */ 145 ierr = MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);CHKERRQ(ierr); 146 147 ierr = ISDestroy(&iscomp);CHKERRQ(ierr); 148 break; 149 case TAO_SUBSET_MATRIXFREE: 150 ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr); 151 ierr = MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);CHKERRQ(ierr); 152 ierr = ISDestroy(&iscomp);CHKERRQ(ierr); 153 break; 154 } 155 PetscFunctionReturn(0); 156 } 157