1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 PETSC_EXTERN PetscErrorCode VecGetRootType_Private(Vec,VecType*); 5 6 typedef struct { 7 Mat A; /* sparse matrix */ 8 Mat U,V; /* dense tall-skinny matrices */ 9 Vec c; /* sequential vector containing the diagonal of C */ 10 Vec work1,work2; /* sequential vectors that hold partial products */ 11 Vec xl,yl; /* auxiliary sequential vectors for matmult operation */ 12 } Mat_LRC; 13 14 static PetscErrorCode MatMult_LRC_kernel(Mat N,Vec x,Vec y,PetscBool transpose) 15 { 16 Mat_LRC *Na = (Mat_LRC*)N->data; 17 PetscErrorCode ierr; 18 PetscMPIInt size; 19 Mat U,V; 20 21 PetscFunctionBegin; 22 U = transpose ? Na->V : Na->U; 23 V = transpose ? Na->U : Na->V; 24 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)N),&size);CHKERRMPI(ierr); 25 if (size == 1) { 26 ierr = MatMultHermitianTranspose(V,x,Na->work1);CHKERRQ(ierr); 27 if (Na->c) { 28 ierr = VecPointwiseMult(Na->work1,Na->c,Na->work1);CHKERRQ(ierr); 29 } 30 if (Na->A) { 31 if (transpose) { 32 ierr = MatMultTranspose(Na->A,x,y);CHKERRQ(ierr); 33 } else { 34 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 35 } 36 ierr = MatMultAdd(U,Na->work1,y,y);CHKERRQ(ierr); 37 } else { 38 ierr = MatMult(U,Na->work1,y);CHKERRQ(ierr); 39 } 40 } else { 41 Mat Uloc,Vloc; 42 Vec yl,xl; 43 const PetscScalar *w1; 44 PetscScalar *w2; 45 PetscInt nwork; 46 PetscMPIInt mpinwork; 47 48 xl = transpose ? Na->yl : Na->xl; 49 yl = transpose ? Na->xl : Na->yl; 50 ierr = VecGetLocalVector(y,yl);CHKERRQ(ierr); 51 ierr = MatDenseGetLocalMatrix(U,&Uloc);CHKERRQ(ierr); 52 ierr = MatDenseGetLocalMatrix(V,&Vloc);CHKERRQ(ierr); 53 54 /* multiply the local part of V with the local part of x */ 55 ierr = VecGetLocalVectorRead(x,xl);CHKERRQ(ierr); 56 ierr = MatMultHermitianTranspose(Vloc,xl,Na->work1);CHKERRQ(ierr); 57 ierr = VecRestoreLocalVectorRead(x,xl);CHKERRQ(ierr); 58 59 /* form the sum of all the local multiplies: this is work2 = V'*x = 60 sum_{all processors} work1 */ 61 ierr = VecGetArrayRead(Na->work1,&w1);CHKERRQ(ierr); 62 ierr = VecGetArrayWrite(Na->work2,&w2);CHKERRQ(ierr); 63 ierr = VecGetLocalSize(Na->work1,&nwork);CHKERRQ(ierr); 64 ierr = PetscMPIIntCast(nwork,&mpinwork);CHKERRQ(ierr); 65 ierr = MPIU_Allreduce(w1,w2,mpinwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRMPI(ierr); 66 ierr = VecRestoreArrayRead(Na->work1,&w1);CHKERRQ(ierr); 67 ierr = VecRestoreArrayWrite(Na->work2,&w2);CHKERRQ(ierr); 68 69 if (Na->c) { /* work2 = C*work2 */ 70 ierr = VecPointwiseMult(Na->work2,Na->c,Na->work2);CHKERRQ(ierr); 71 } 72 73 if (Na->A) { 74 /* form y = A*x or A^t*x */ 75 if (transpose) { 76 ierr = MatMultTranspose(Na->A,x,y);CHKERRQ(ierr); 77 } else { 78 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 79 } 80 /* multiply-add y = y + U*work2 */ 81 ierr = MatMultAdd(Uloc,Na->work2,yl,yl);CHKERRQ(ierr); 82 } else { 83 /* multiply y = U*work2 */ 84 ierr = MatMult(Uloc,Na->work2,yl);CHKERRQ(ierr); 85 } 86 87 ierr = VecRestoreLocalVector(y,yl);CHKERRQ(ierr); 88 } 89 PetscFunctionReturn(0); 90 } 91 92 static PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y) 93 { 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = MatMult_LRC_kernel(N,x,y,PETSC_FALSE);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 static PetscErrorCode MatMultTranspose_LRC(Mat N,Vec x,Vec y) 102 { 103 PetscErrorCode ierr; 104 105 PetscFunctionBegin; 106 ierr = MatMult_LRC_kernel(N,x,y,PETSC_TRUE);CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode MatDestroy_LRC(Mat N) 111 { 112 Mat_LRC *Na = (Mat_LRC*)N->data; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = MatDestroy(&Na->A);CHKERRQ(ierr); 117 ierr = MatDestroy(&Na->U);CHKERRQ(ierr); 118 ierr = MatDestroy(&Na->V);CHKERRQ(ierr); 119 ierr = VecDestroy(&Na->c);CHKERRQ(ierr); 120 ierr = VecDestroy(&Na->work1);CHKERRQ(ierr); 121 ierr = VecDestroy(&Na->work2);CHKERRQ(ierr); 122 ierr = VecDestroy(&Na->xl);CHKERRQ(ierr); 123 ierr = VecDestroy(&Na->yl);CHKERRQ(ierr); 124 ierr = PetscFree(N->data);CHKERRQ(ierr); 125 ierr = PetscObjectComposeFunction((PetscObject)N,"MatLRCGetMats_C",NULL);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 static PetscErrorCode MatLRCGetMats_LRC(Mat N,Mat *A,Mat *U,Vec *c,Mat *V) 130 { 131 Mat_LRC *Na = (Mat_LRC*)N->data; 132 133 PetscFunctionBegin; 134 if (A) *A = Na->A; 135 if (U) *U = Na->U; 136 if (c) *c = Na->c; 137 if (V) *V = Na->V; 138 PetscFunctionReturn(0); 139 } 140 141 /*@ 142 MatLRCGetMats - Returns the constituents of an LRC matrix 143 144 Collective on Mat 145 146 Input Parameter: 147 . N - matrix of type LRC 148 149 Output Parameters: 150 + A - the (sparse) matrix 151 . U - first dense rectangular (tall and skinny) matrix 152 . c - a sequential vector containing the diagonal of C 153 - V - second dense rectangular (tall and skinny) matrix 154 155 Note: 156 The returned matrices need not be destroyed by the caller. 157 158 Level: intermediate 159 160 .seealso: MatCreateLRC() 161 @*/ 162 PetscErrorCode MatLRCGetMats(Mat N,Mat *A,Mat *U,Vec *c,Mat *V) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 ierr = PetscUseMethod(N,"MatLRCGetMats_C",(Mat,Mat*,Mat*,Vec*,Mat*),(N,A,U,c,V));CHKERRQ(ierr); 168 PetscFunctionReturn(0); 169 } 170 171 /*@ 172 MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V' 173 174 Collective on Mat 175 176 Input Parameters: 177 + A - the (sparse) matrix (can be NULL) 178 . U, V - two dense rectangular (tall and skinny) matrices 179 - c - a vector containing the diagonal of C (can be NULL) 180 181 Output Parameter: 182 . N - the matrix that represents A + U*C*V' 183 184 Notes: 185 The matrix A + U*C*V' is not formed! Rather the new matrix 186 object performs the matrix-vector product by first multiplying by 187 A and then adding the other term. 188 189 C is a diagonal matrix (represented as a vector) of order k, 190 where k is the number of columns of both U and V. 191 192 If A is NULL then the new object behaves like a low-rank matrix U*C*V'. 193 194 Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*C*U'. 195 196 If c is NULL then the low-rank correction is just U*V'. 197 If a sequential c vector is used for a parallel matrix, 198 PETSc assumes that the values of the vector are consistently set across processors. 199 200 Level: intermediate 201 202 .seealso: MatLRCGetMats() 203 @*/ 204 PetscErrorCode MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat *N) 205 { 206 PetscErrorCode ierr; 207 PetscBool match; 208 PetscInt m,n,k,m1,n1,k1; 209 Mat_LRC *Na; 210 Mat Uloc; 211 PetscMPIInt size, csize = 0; 212 213 PetscFunctionBegin; 214 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,1); 215 PetscValidHeaderSpecific(U,MAT_CLASSID,2); 216 if (c) PetscValidHeaderSpecific(c,VEC_CLASSID,3); 217 if (V) { 218 PetscValidHeaderSpecific(V,MAT_CLASSID,4); 219 PetscCheckSameComm(U,2,V,4); 220 } 221 if (A) PetscCheckSameComm(A,1,U,2); 222 223 if (!V) V = U; 224 ierr = PetscObjectBaseTypeCompareAny((PetscObject)U,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 225 PetscCheckFalse(!match,PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix U must be of type dense, found %s",((PetscObject)U)->type_name); 226 ierr = PetscObjectBaseTypeCompareAny((PetscObject)V,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 227 PetscCheckFalse(!match,PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix V must be of type dense, found %s",((PetscObject)V)->type_name); 228 ierr = PetscStrcmp(U->defaultvectype,V->defaultvectype,&match);CHKERRQ(ierr); 229 PetscCheckFalse(!match,PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_WRONG,"Matrix U and V must have the same VecType %s != %s",U->defaultvectype,V->defaultvectype); 230 if (A) { 231 ierr = PetscStrcmp(A->defaultvectype,U->defaultvectype,&match);CHKERRQ(ierr); 232 PetscCheckFalse(!match,PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_WRONG,"Matrix A and U must have the same VecType %s != %s",A->defaultvectype,U->defaultvectype); 233 } 234 235 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)U),&size);CHKERRMPI(ierr); 236 ierr = MatGetSize(U,NULL,&k);CHKERRQ(ierr); 237 ierr = MatGetSize(V,NULL,&k1);CHKERRQ(ierr); 238 PetscCheckFalse(k != k1,PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"U and V have different number of columns (%" PetscInt_FMT " vs %" PetscInt_FMT ")",k,k1); 239 ierr = MatGetLocalSize(U,&m,NULL);CHKERRQ(ierr); 240 ierr = MatGetLocalSize(V,&n,NULL);CHKERRQ(ierr); 241 if (A) { 242 ierr = MatGetLocalSize(A,&m1,&n1);CHKERRQ(ierr); 243 PetscCheckFalse(m != m1,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local dimensions of U %" PetscInt_FMT " and A %" PetscInt_FMT " do not match",m,m1); 244 PetscCheckFalse(n != n1,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local dimensions of V %" PetscInt_FMT " and A %" PetscInt_FMT " do not match",n,n1); 245 } 246 if (c) { 247 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)c),&csize);CHKERRMPI(ierr); 248 ierr = VecGetSize(c,&k1);CHKERRQ(ierr); 249 PetscCheckFalse(k != k1,PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"The length of c %" PetscInt_FMT " does not match the number of columns of U and V (%" PetscInt_FMT ")",k1,k); 250 PetscCheckFalse(csize != 1 && csize != size, PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"U and c must have the same communicator size %d != %d",size,csize); 251 } 252 253 ierr = MatCreate(PetscObjectComm((PetscObject)U),N);CHKERRQ(ierr); 254 ierr = MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 255 ierr = MatSetVecType(*N,U->defaultvectype);CHKERRQ(ierr); 256 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATLRC);CHKERRQ(ierr); 257 /* Flag matrix as symmetric if A is symmetric and U == V */ 258 ierr = MatSetOption(*N,MAT_SYMMETRIC,(PetscBool)((A ? A->symmetric : PETSC_TRUE) && U == V));CHKERRQ(ierr); 259 260 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 261 (*N)->data = (void*)Na; 262 Na->A = A; 263 Na->U = U; 264 Na->c = c; 265 Na->V = V; 266 267 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 268 ierr = PetscObjectReference((PetscObject)Na->U);CHKERRQ(ierr); 269 ierr = PetscObjectReference((PetscObject)Na->V);CHKERRQ(ierr); 270 ierr = PetscObjectReference((PetscObject)c);CHKERRQ(ierr); 271 272 ierr = MatDenseGetLocalMatrix(Na->U,&Uloc);CHKERRQ(ierr); 273 ierr = MatCreateVecs(Uloc,&Na->work1,NULL);CHKERRQ(ierr); 274 if (size != 1) { 275 Mat Vloc; 276 277 if (Na->c && csize != 1) { /* scatter parallel vector to sequential */ 278 VecScatter sct; 279 280 ierr = VecScatterCreateToAll(Na->c,&sct,&c);CHKERRQ(ierr); 281 ierr = VecScatterBegin(sct,Na->c,c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 282 ierr = VecScatterEnd(sct,Na->c,c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 283 ierr = VecScatterDestroy(&sct);CHKERRQ(ierr); 284 ierr = VecDestroy(&Na->c);CHKERRQ(ierr); 285 ierr = PetscLogObjectParent((PetscObject)*N,(PetscObject)c);CHKERRQ(ierr); 286 Na->c = c; 287 } 288 ierr = MatDenseGetLocalMatrix(Na->V,&Vloc);CHKERRQ(ierr); 289 ierr = VecDuplicate(Na->work1,&Na->work2);CHKERRQ(ierr); 290 ierr = MatCreateVecs(Vloc,NULL,&Na->xl);CHKERRQ(ierr); 291 ierr = MatCreateVecs(Uloc,NULL,&Na->yl);CHKERRQ(ierr); 292 } 293 ierr = PetscLogObjectParent((PetscObject)*N,(PetscObject)Na->work1);CHKERRQ(ierr); 294 ierr = PetscLogObjectParent((PetscObject)*N,(PetscObject)Na->work1);CHKERRQ(ierr); 295 ierr = PetscLogObjectParent((PetscObject)*N,(PetscObject)Na->xl);CHKERRQ(ierr); 296 ierr = PetscLogObjectParent((PetscObject)*N,(PetscObject)Na->yl);CHKERRQ(ierr); 297 298 /* Internally create a scaling vector if roottypes do not match */ 299 if (Na->c) { 300 VecType rt1,rt2; 301 302 ierr = VecGetRootType_Private(Na->work1,&rt1);CHKERRQ(ierr); 303 ierr = VecGetRootType_Private(Na->c,&rt2);CHKERRQ(ierr); 304 ierr = PetscStrcmp(rt1,rt2,&match);CHKERRQ(ierr); 305 if (!match) { 306 ierr = VecDuplicate(Na->c,&c);CHKERRQ(ierr); 307 ierr = VecCopy(Na->c,c);CHKERRQ(ierr); 308 ierr = VecDestroy(&Na->c);CHKERRQ(ierr); 309 ierr = PetscLogObjectParent((PetscObject)*N,(PetscObject)c);CHKERRQ(ierr); 310 Na->c = c; 311 } 312 } 313 314 (*N)->ops->destroy = MatDestroy_LRC; 315 (*N)->ops->mult = MatMult_LRC; 316 (*N)->ops->multtranspose = MatMultTranspose_LRC; 317 318 (*N)->assembled = PETSC_TRUE; 319 (*N)->preallocated = PETSC_TRUE; 320 321 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatLRCGetMats_C",MatLRCGetMats_LRC);CHKERRQ(ierr); 322 ierr = MatSetUp(*N);CHKERRQ(ierr); 323 PetscFunctionReturn(0); 324 } 325