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