1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 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; 12 13 14 PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y) 15 { 16 Mat_LRC *Na = (Mat_LRC*)N->data; 17 PetscErrorCode ierr; 18 PetscScalar *w1,*w2; 19 const PetscScalar *a; 20 21 PetscFunctionBegin; 22 ierr = VecGetArrayRead(x,&a);CHKERRQ(ierr); 23 ierr = VecPlaceArray(Na->xl,a);CHKERRQ(ierr); 24 ierr = VecGetLocalVector(y,Na->yl);CHKERRQ(ierr); 25 26 /* multiply the local part of V with the local part of x */ 27 #if defined(PETSC_USE_COMPLEX) 28 ierr = MatMultHermitianTranspose(Na->V,Na->xl,Na->work1);CHKERRQ(ierr); 29 #else 30 ierr = MatMultTranspose(Na->V,Na->xl,Na->work1);CHKERRQ(ierr); 31 #endif 32 33 /* form the sum of all the local multiplies: this is work2 = V'*x = 34 sum_{all processors} work1 */ 35 ierr = VecGetArray(Na->work1,&w1);CHKERRQ(ierr); 36 ierr = VecGetArray(Na->work2,&w2);CHKERRQ(ierr); 37 ierr = MPIU_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRQ(ierr); 38 ierr = VecRestoreArray(Na->work1,&w1);CHKERRQ(ierr); 39 ierr = VecRestoreArray(Na->work2,&w2);CHKERRQ(ierr); 40 41 if (Na->c) { /* work2 = C*work2 */ 42 ierr = VecPointwiseMult(Na->work2,Na->c,Na->work2);CHKERRQ(ierr); 43 } 44 45 if (Na->A) { 46 /* form y = A*x */ 47 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 48 /* multiply-add y = y + U*work2 */ 49 ierr = MatMultAdd(Na->U,Na->work2,Na->yl,Na->yl);CHKERRQ(ierr); 50 } else { 51 /* multiply y = U*work2 */ 52 ierr = MatMult(Na->U,Na->work2,Na->yl);CHKERRQ(ierr); 53 } 54 55 ierr = VecRestoreArrayRead(x,&a);CHKERRQ(ierr); 56 ierr = VecResetArray(Na->xl);CHKERRQ(ierr); 57 ierr = VecRestoreLocalVector(y,Na->yl);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 PetscErrorCode MatDestroy_LRC(Mat N) 62 { 63 Mat_LRC *Na = (Mat_LRC*)N->data; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 ierr = MatDestroy(&Na->A);CHKERRQ(ierr); 68 ierr = MatDestroy(&Na->U);CHKERRQ(ierr); 69 ierr = MatDestroy(&Na->V);CHKERRQ(ierr); 70 ierr = VecDestroy(&Na->c);CHKERRQ(ierr); 71 ierr = VecDestroy(&Na->work1);CHKERRQ(ierr); 72 ierr = VecDestroy(&Na->work2);CHKERRQ(ierr); 73 ierr = VecDestroy(&Na->xl);CHKERRQ(ierr); 74 ierr = VecDestroy(&Na->yl);CHKERRQ(ierr); 75 ierr = PetscFree(N->data);CHKERRQ(ierr); 76 ierr = PetscObjectComposeFunction((PetscObject)N,"MatLRCGetMats_C",0);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 PetscErrorCode MatLRCGetMats_LRC(Mat N,Mat *A,Mat *U,Vec *c,Mat *V) 81 { 82 Mat_LRC *Na = (Mat_LRC*)N->data; 83 84 PetscFunctionBegin; 85 if (A) *A = Na->A; 86 if (U) *U = Na->U; 87 if (c) *c = Na->c; 88 if (V) *V = Na->V; 89 PetscFunctionReturn(0); 90 } 91 92 /*@ 93 MatLRCGetMats - Returns the constituents of an LRC matrix 94 95 Collective on Mat 96 97 Input Parameter: 98 . N - matrix of type LRC 99 100 Output Parameters: 101 + A - the (sparse) matrix 102 . U, V - two dense rectangular (tall and skinny) matrices 103 - c - a sequential vector containing the diagonal of C 104 105 Note: 106 The returned matrices need not be destroyed by the caller. 107 108 Level: intermediate 109 110 .seealso: MatCreateLRC() 111 @*/ 112 PetscErrorCode MatLRCGetMats(Mat N,Mat *A,Mat *U,Vec *c,Mat *V) 113 { 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 ierr = PetscUseMethod(N,"MatLRCGetMats_C",(Mat,Mat*,Mat*,Vec*,Mat*),(N,A,U,c,V));CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 /*@ 122 MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V' 123 124 Collective on Mat 125 126 Input Parameters: 127 + A - the (sparse) matrix (can be NULL) 128 . U, V - two dense rectangular (tall and skinny) matrices 129 - c - a sequential vector containing the diagonal of C (can be NULL) 130 131 Output Parameter: 132 . N - the matrix that represents A + U*C*V' 133 134 Notes: 135 The matrix A + U*C*V' is not formed! Rather the new matrix 136 object performs the matrix-vector product by first multiplying by 137 A and then adding the other term. 138 139 C is a diagonal matrix (represented as a vector) of order k, 140 where k is the number of columns of both U and V. 141 142 If A is NULL then the new object behaves like a low-rank matrix U*C*V'. 143 144 Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*C*U'. 145 146 If c is NULL then the low-rank correction is just U*V'. 147 148 Level: intermediate 149 150 .seealso: MatLRCGetMats() 151 @*/ 152 PetscErrorCode MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat *N) 153 { 154 PetscErrorCode ierr; 155 PetscBool match; 156 PetscInt m,n,k,m1,n1,k1; 157 Mat_LRC *Na; 158 159 PetscFunctionBegin; 160 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,1); 161 PetscValidHeaderSpecific(U,MAT_CLASSID,2); 162 if (c) PetscValidHeaderSpecific(c,VEC_CLASSID,3); 163 if (V) PetscValidHeaderSpecific(V,MAT_CLASSID,4); 164 else V=U; 165 if (A) PetscCheckSameComm(A,1,U,2); 166 PetscCheckSameComm(U,2,V,4); 167 168 ierr = PetscObjectTypeCompareAny((PetscObject)U,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 169 if (!match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix U must be of type dense"); 170 if (V) { 171 ierr = PetscObjectTypeCompareAny((PetscObject)V,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 172 if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Matrix V must be of type dense"); 173 } 174 175 ierr = MatGetSize(U,NULL,&k);CHKERRQ(ierr); 176 ierr = MatGetSize(V,NULL,&k1);CHKERRQ(ierr); 177 if (k!=k1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"U and V have different number of columns (%D vs %D)",k,k1); 178 ierr = MatGetLocalSize(U,&m,NULL);CHKERRQ(ierr); 179 ierr = MatGetLocalSize(V,&n,NULL);CHKERRQ(ierr); 180 if (A) { 181 ierr = MatGetLocalSize(A,&m1,&n1);CHKERRQ(ierr); 182 if (m!=m1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"Local dimensions of U %D and A %D do not match",m,m1); 183 if (n!=n1) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Local dimensions of V %D and A %D do not match",n,n1); 184 } 185 if (c) { 186 ierr = VecGetSize(c,&k1);CHKERRQ(ierr); 187 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); 188 ierr = VecGetLocalSize(c,&k1);CHKERRQ(ierr); 189 if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"c must be a sequential vector"); 190 } 191 192 ierr = MatCreate(PetscObjectComm((PetscObject)U),N);CHKERRQ(ierr); 193 ierr = MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 194 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATLRC);CHKERRQ(ierr); 195 196 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 197 (*N)->data = (void*)Na; 198 Na->A = A; 199 Na->c = c; 200 201 ierr = MatDenseGetLocalMatrix(U,&Na->U);CHKERRQ(ierr); 202 ierr = MatDenseGetLocalMatrix(V,&Na->V);CHKERRQ(ierr); 203 if (A) { ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); } 204 ierr = PetscObjectReference((PetscObject)Na->U);CHKERRQ(ierr); 205 ierr = PetscObjectReference((PetscObject)Na->V);CHKERRQ(ierr); 206 if (c) { ierr = PetscObjectReference((PetscObject)c);CHKERRQ(ierr); } 207 208 ierr = VecCreateSeq(PETSC_COMM_SELF,U->cmap->N,&Na->work1);CHKERRQ(ierr); 209 ierr = VecDuplicate(Na->work1,&Na->work2);CHKERRQ(ierr); 210 ierr = PetscMPIIntCast(U->cmap->N,&Na->nwork);CHKERRQ(ierr); 211 212 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,V->rmap->n,NULL,&Na->xl);CHKERRQ(ierr); 213 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,U->rmap->n,NULL,&Na->yl);CHKERRQ(ierr); 214 215 (*N)->ops->destroy = MatDestroy_LRC; 216 (*N)->ops->mult = MatMult_LRC; 217 (*N)->assembled = PETSC_TRUE; 218 (*N)->preallocated = PETSC_TRUE; 219 (*N)->cmap->N = V->rmap->N; 220 (*N)->rmap->N = U->rmap->N; 221 (*N)->cmap->n = V->rmap->n; 222 (*N)->rmap->n = U->rmap->n; 223 224 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatLRCGetMats_C",MatLRCGetMats_LRC);CHKERRQ(ierr); 225 PetscFunctionReturn(0); 226 } 227 228