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>

  8: /*@C
  9:   TaoVecGetSubVec - Gets a subvector using the IS

 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)

 17:   Output Parameters:
 18: . vreduced - the subvector

 20:   Notes:
 21:   maskvalue should usually be 0.0, unless a pointwise divide will be used.

 23: @*/
 24: PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced)
 25: {
 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;


 40:   VecGetSize(vfull, &nfull);
 41:   ISGetSize(is, &nreduced);

 43:   if (nreduced == nfull) {
 44:     VecDestroy(vreduced);
 45:     VecDuplicate(vfull,vreduced);
 46:     VecCopy(vfull,*vreduced);
 47:   } else {
 48:     switch (reduced_type) {
 49:     case TAO_SUBSET_SUBVEC:
 50:       VecGetType(vfull,&vtype);
 51:       VecGetOwnershipRange(vfull,&flow,&fhigh);
 52:       ISGetLocalSize(is,&nreduced_local);
 53:       PetscObjectGetComm((PetscObject)vfull,&comm);
 54:       if (*vreduced) {
 55:         VecDestroy(vreduced);
 56:       }
 57:       VecCreate(comm,vreduced);
 58:       VecSetType(*vreduced,vtype);

 60:       VecSetSizes(*vreduced,nreduced_local,nreduced);
 61:       VecGetOwnershipRange(*vreduced,&rlow,&rhigh);
 62:       ISCreateStride(comm,nreduced_local,rlow,1,&ident);
 63:       VecScatterCreate(vfull,is,*vreduced,ident,&scatter);
 64:       VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
 65:       VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
 66:       VecScatterDestroy(&scatter);
 67:       ISDestroy(&ident);
 68:       break;

 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:         VecDuplicate(vfull,vreduced);
 76:       }

 78:       VecSet(*vreduced,maskvalue);
 79:       ISGetLocalSize(is,&nlocal);
 80:       VecGetOwnershipRange(vfull,&flow,&fhigh);
 81:       VecGetArray(vfull,&fv);
 82:       VecGetArray(*vreduced,&rv);
 83:       ISGetIndices(is,&s);
 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:       ISRestoreIndices(is,&s);
 89:       VecRestoreArray(vfull,&fv);
 90:       VecRestoreArray(*vreduced,&rv);
 91:       break;
 92:     }
 93:   }
 94:   return(0);
 95: }

 99: /*@C
100:   TaoMatGetSubMat - Gets a submatrix using the IS

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)

109:   Output Parameters:
110: . Msub - the submatrix
111: @*/
112: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
113: {
115:   IS             iscomp;
116:   PetscBool      flg = PETSC_FALSE;

121:   MatDestroy(Msub);
122:   switch (subset_type) {
123:   case TAO_SUBSET_SUBVEC:
124:     MatGetSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
125:     break;

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:     PetscOptionsBegin(PetscObjectComm((PetscObject)M),NULL,NULL,NULL);
133:     PetscOptionsBool("-different_submatrix","use separate hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);
134:     PetscOptionsEnd();
135:     if (flg) {
136:       MatDuplicate(M, MAT_COPY_VALUES, Msub);
137:     } else {
138:       /* Act on hessian directly (default) */
139:       PetscObjectReference((PetscObject)M);
140:       *Msub = M;
141:     }
142:     /* Save the diagonal to temporary vector */
143:     MatGetDiagonal(*Msub,v1);

145:     /* Zero out rows and columns */
146:     ISComplementVec(is,v1,&iscomp);

148:     /* Use v1 instead of 0 here because of PETSc bug */
149:     MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);

151:     ISDestroy(&iscomp);
152:     break;
153:   case TAO_SUBSET_MATRIXFREE:
154:     ISComplementVec(is,v1,&iscomp);
155:     MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);
156:     ISDestroy(&iscomp);
157:     break;
158:   }
159:   return(0);
160: }