Actual source code: lrc.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  2:  #include <petsc/private/matimpl.h>

  4: typedef struct {
  5:   Mat         A;           /* sparse matrix */
  6:   Mat         U,V;         /* dense tall-skinny matrices */
  7:   Vec         c;           /* sequential vector containing the diagonal of C */
  8:   Vec         work1,work2; /* sequential vectors that hold partial products */
  9:   PetscMPIInt nwork;       /* length of work vectors */
 10:   Vec         xl,yl;       /* auxiliary sequential vectors for matmult operation */
 11: } Mat_LRC;


 14: PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y)
 15: {
 16:   Mat_LRC           *Na = (Mat_LRC*)N->data;
 17:   Mat               Uloc,Vloc;
 18:   PetscErrorCode    ierr;
 19:   PetscScalar       *w1,*w2;
 20:   const PetscScalar *a;

 23:   VecGetArrayRead(x,&a);
 24:   VecPlaceArray(Na->xl,a);
 25:   VecGetLocalVector(y,Na->yl);
 26:   MatDenseGetLocalMatrix(Na->U,&Uloc);
 27:   MatDenseGetLocalMatrix(Na->V,&Vloc);

 29:   /* multiply the local part of V with the local part of x */
 30: #if defined(PETSC_USE_COMPLEX)
 31:   MatMultHermitianTranspose(Vloc,Na->xl,Na->work1);
 32: #else
 33:   MatMultTranspose(Vloc,Na->xl,Na->work1);
 34: #endif

 36:   /* form the sum of all the local multiplies: this is work2 = V'*x =
 37:      sum_{all processors} work1 */
 38:   VecGetArray(Na->work1,&w1);
 39:   VecGetArray(Na->work2,&w2);
 40:   MPIU_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
 41:   VecRestoreArray(Na->work1,&w1);
 42:   VecRestoreArray(Na->work2,&w2);

 44:   if (Na->c) {  /* work2 = C*work2 */
 45:     VecPointwiseMult(Na->work2,Na->c,Na->work2);
 46:   }

 48:   if (Na->A) {
 49:     /* form y = A*x */
 50:     MatMult(Na->A,x,y);
 51:     /* multiply-add y = y + U*work2 */
 52:     MatMultAdd(Uloc,Na->work2,Na->yl,Na->yl);
 53:   } else {
 54:     /* multiply y = U*work2 */
 55:     MatMult(Uloc,Na->work2,Na->yl);
 56:   }

 58:   VecRestoreArrayRead(x,&a);
 59:   VecResetArray(Na->xl);
 60:   VecRestoreLocalVector(y,Na->yl);
 61:   return(0);
 62: }

 64: PetscErrorCode MatDestroy_LRC(Mat N)
 65: {
 66:   Mat_LRC        *Na = (Mat_LRC*)N->data;

 70:   MatDestroy(&Na->A);
 71:   MatDestroy(&Na->U);
 72:   MatDestroy(&Na->V);
 73:   VecDestroy(&Na->c);
 74:   VecDestroy(&Na->work1);
 75:   VecDestroy(&Na->work2);
 76:   VecDestroy(&Na->xl);
 77:   VecDestroy(&Na->yl);
 78:   PetscFree(N->data);
 79:   PetscObjectComposeFunction((PetscObject)N,"MatLRCGetMats_C",0);
 80:   return(0);
 81: }

 83: PetscErrorCode MatLRCGetMats_LRC(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
 84: {
 85:   Mat_LRC *Na = (Mat_LRC*)N->data;

 88:   if (A) *A = Na->A;
 89:   if (U) *U = Na->U;
 90:   if (c) *c = Na->c;
 91:   if (V) *V = Na->V;
 92:   return(0);
 93: }

 95: /*@
 96:    MatLRCGetMats - Returns the constituents of an LRC matrix

 98:    Collective on Mat

100:    Input Parameter:
101: .  N    - matrix of type LRC

103:    Output Parameters:
104: +  A    - the (sparse) matrix
105: .  U, V - two dense rectangular (tall and skinny) matrices
106: -  c    - a sequential vector containing the diagonal of C

108:    Note:
109:    The returned matrices need not be destroyed by the caller.

111:    Level: intermediate

113: .seealso: MatCreateLRC()
114: @*/
115: PetscErrorCode  MatLRCGetMats(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
116: {

120:   PetscUseMethod(N,"MatLRCGetMats_C",(Mat,Mat*,Mat*,Vec*,Mat*),(N,A,U,c,V));
121:   return(0);
122: }

124: /*@
125:    MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V'

127:    Collective on Mat

129:    Input Parameters:
130: +  A    - the (sparse) matrix (can be NULL)
131: .  U, V - two dense rectangular (tall and skinny) matrices
132: -  c    - a sequential vector containing the diagonal of C (can be NULL) 

134:    Output Parameter:
135: .  N    - the matrix that represents A + U*C*V'

137:    Notes:
138:    The matrix A + U*C*V' is not formed! Rather the new matrix
139:    object performs the matrix-vector product by first multiplying by
140:    A and then adding the other term.

142:    C is a diagonal matrix (represented as a vector) of order k,
143:    where k is the number of columns of both U and V.

145:    If A is NULL then the new object behaves like a low-rank matrix U*C*V'.

147:    Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*C*U'.

149:    If c is NULL then the low-rank correction is just U*V'.

151:    Level: intermediate

153: .seealso: MatLRCGetMats()
154: @*/
155: PetscErrorCode MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat *N)
156: {
158:   PetscBool      match;
159:   PetscInt       m,n,k,m1,n1,k1;
160:   Mat_LRC        *Na;

167:   else V=U;

171:   PetscObjectTypeCompareAny((PetscObject)U,&match,MATSEQDENSE,MATMPIDENSE,"");
172:   if (!match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix U must be of type dense");
173:   if (V) {
174:     PetscObjectTypeCompareAny((PetscObject)V,&match,MATSEQDENSE,MATMPIDENSE,"");
175:     if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Matrix V must be of type dense");
176:   }

178:   MatGetSize(U,NULL,&k);
179:   MatGetSize(V,NULL,&k1);
180:   if (k!=k1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"U and V have different number of columns (%D vs %D)",k,k1);
181:   MatGetLocalSize(U,&m,NULL);
182:   MatGetLocalSize(V,&n,NULL);
183:   if (A) {
184:     MatGetLocalSize(A,&m1,&n1);
185:     if (m!=m1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"Local dimensions of U %D and A %D do not match",m,m1);
186:     if (n!=n1) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Local dimensions of V %D and A %D do not match",n,n1);
187:   }
188:   if (c) {
189:     VecGetSize(c,&k1);
190:     if (k!=k1) SETERRQ2(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"The length of c %D does not match the number of columns of U and V (%D)",k1,k);
191:     VecGetLocalSize(c,&k1);
192:     if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"c must be a sequential vector");
193:   }

195:   MatCreate(PetscObjectComm((PetscObject)U),N);
196:   MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);
197:   PetscObjectChangeTypeName((PetscObject)*N,MATLRC);

199:   PetscNewLog(*N,&Na);
200:   (*N)->data = (void*)Na;
201:   Na->A      = A;
202:   Na->U      = U;
203:   Na->c      = c;
204:   Na->V      = V;

206:   if (A) { PetscObjectReference((PetscObject)A); }
207:   PetscObjectReference((PetscObject)Na->U);
208:   PetscObjectReference((PetscObject)Na->V);
209:   if (c) { PetscObjectReference((PetscObject)c); }

211:   VecCreateSeq(PETSC_COMM_SELF,U->cmap->N,&Na->work1);
212:   VecDuplicate(Na->work1,&Na->work2);
213:   PetscMPIIntCast(U->cmap->N,&Na->nwork);

215:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,V->rmap->n,NULL,&Na->xl);
216:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,U->rmap->n,NULL,&Na->yl);

218:   (*N)->ops->destroy = MatDestroy_LRC;
219:   (*N)->ops->mult    = MatMult_LRC;
220:   (*N)->assembled    = PETSC_TRUE;
221:   (*N)->preallocated = PETSC_TRUE;
222:   (*N)->cmap->N      = V->rmap->N;
223:   (*N)->rmap->N      = U->rmap->N;
224:   (*N)->cmap->n      = V->rmap->n;
225:   (*N)->rmap->n      = U->rmap->n;

227:   PetscObjectComposeFunction((PetscObject)(*N),"MatLRCGetMats_C",MatLRCGetMats_LRC);
228:   return(0);
229: }