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