1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; 6 Mat D; /* local submatrix for diagonal part */ 7 Vec w, left, right, leftwork, rightwork; 8 PetscScalar scale; 9 } Mat_Normal; 10 11 PetscErrorCode MatScale_NormalHermitian(Mat inA, PetscScalar scale) { 12 Mat_Normal *a = (Mat_Normal *)inA->data; 13 14 PetscFunctionBegin; 15 a->scale *= scale; 16 PetscFunctionReturn(0); 17 } 18 19 PetscErrorCode MatDiagonalScale_NormalHermitian(Mat inA, Vec left, Vec right) { 20 Mat_Normal *a = (Mat_Normal *)inA->data; 21 22 PetscFunctionBegin; 23 if (left) { 24 if (!a->left) { 25 PetscCall(VecDuplicate(left, &a->left)); 26 PetscCall(VecCopy(left, a->left)); 27 } else { 28 PetscCall(VecPointwiseMult(a->left, left, a->left)); 29 } 30 } 31 if (right) { 32 if (!a->right) { 33 PetscCall(VecDuplicate(right, &a->right)); 34 PetscCall(VecCopy(right, a->right)); 35 } else { 36 PetscCall(VecPointwiseMult(a->right, right, a->right)); 37 } 38 } 39 PetscFunctionReturn(0); 40 } 41 42 PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) { 43 Mat_Normal *a = (Mat_Normal *)mat->data; 44 Mat B = a->A, *suba; 45 IS *row; 46 PetscInt M; 47 48 PetscFunctionBegin; 49 PetscCheck(!a->left && !a->right && irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented"); 50 if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat)); 51 PetscCall(MatGetSize(B, &M, NULL)); 52 PetscCall(PetscMalloc1(n, &row)); 53 PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0])); 54 PetscCall(ISSetIdentity(row[0])); 55 for (M = 1; M < n; ++M) row[M] = row[0]; 56 PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba)); 57 for (M = 0; M < n; ++M) { 58 PetscCall(MatCreateNormalHermitian(suba[M], *submat + M)); 59 ((Mat_Normal *)(*submat)[M]->data)->scale = a->scale; 60 } 61 PetscCall(ISDestroy(&row[0])); 62 PetscCall(PetscFree(row)); 63 PetscCall(MatDestroySubMatrices(n, &suba)); 64 PetscFunctionReturn(0); 65 } 66 67 PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B) { 68 Mat_Normal *a = (Mat_Normal *)A->data; 69 Mat C, Aa = a->A; 70 IS row; 71 72 PetscFunctionBegin; 73 PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same"); 74 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row)); 75 PetscCall(ISSetIdentity(row)); 76 PetscCall(MatPermute(Aa, row, colp, &C)); 77 PetscCall(ISDestroy(&row)); 78 PetscCall(MatCreateNormalHermitian(C, B)); 79 PetscCall(MatDestroy(&C)); 80 PetscFunctionReturn(0); 81 } 82 83 PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B) { 84 Mat_Normal *a = (Mat_Normal *)A->data; 85 Mat C; 86 87 PetscFunctionBegin; 88 PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented"); 89 PetscCall(MatDuplicate(a->A, op, &C)); 90 PetscCall(MatCreateNormalHermitian(C, B)); 91 PetscCall(MatDestroy(&C)); 92 if (op == MAT_COPY_VALUES) ((Mat_Normal *)(*B)->data)->scale = a->scale; 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str) { 97 Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data; 98 99 PetscFunctionBegin; 100 PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented"); 101 PetscCall(MatCopy(a->A, b->A, str)); 102 b->scale = a->scale; 103 PetscCall(VecDestroy(&b->left)); 104 PetscCall(VecDestroy(&b->right)); 105 PetscCall(VecDestroy(&b->leftwork)); 106 PetscCall(VecDestroy(&b->rightwork)); 107 PetscFunctionReturn(0); 108 } 109 110 PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y) { 111 Mat_Normal *Na = (Mat_Normal *)N->data; 112 Vec in; 113 114 PetscFunctionBegin; 115 in = x; 116 if (Na->right) { 117 if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork)); 118 PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in)); 119 in = Na->rightwork; 120 } 121 PetscCall(MatMult(Na->A, in, Na->w)); 122 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y)); 123 if (Na->left) PetscCall(VecPointwiseMult(y, Na->left, y)); 124 PetscCall(VecScale(y, Na->scale)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode MatMultHermitianAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3) { 129 Mat_Normal *Na = (Mat_Normal *)N->data; 130 Vec in; 131 132 PetscFunctionBegin; 133 in = v1; 134 if (Na->right) { 135 if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork)); 136 PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in)); 137 in = Na->rightwork; 138 } 139 PetscCall(MatMult(Na->A, in, Na->w)); 140 PetscCall(VecScale(Na->w, Na->scale)); 141 if (Na->left) { 142 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3)); 143 PetscCall(VecPointwiseMult(v3, Na->left, v3)); 144 PetscCall(VecAXPY(v3, 1.0, v2)); 145 } else { 146 PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3)); 147 } 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MatMultHermitianTranspose_Normal(Mat N, Vec x, Vec y) { 152 Mat_Normal *Na = (Mat_Normal *)N->data; 153 Vec in; 154 155 PetscFunctionBegin; 156 in = x; 157 if (Na->left) { 158 if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork)); 159 PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in)); 160 in = Na->leftwork; 161 } 162 PetscCall(MatMult(Na->A, in, Na->w)); 163 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y)); 164 if (Na->right) PetscCall(VecPointwiseMult(y, Na->right, y)); 165 PetscCall(VecScale(y, Na->scale)); 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3) { 170 Mat_Normal *Na = (Mat_Normal *)N->data; 171 Vec in; 172 173 PetscFunctionBegin; 174 in = v1; 175 if (Na->left) { 176 if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork)); 177 PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in)); 178 in = Na->leftwork; 179 } 180 PetscCall(MatMult(Na->A, in, Na->w)); 181 PetscCall(VecScale(Na->w, Na->scale)); 182 if (Na->right) { 183 PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3)); 184 PetscCall(VecPointwiseMult(v3, Na->right, v3)); 185 PetscCall(VecAXPY(v3, 1.0, v2)); 186 } else { 187 PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3)); 188 } 189 PetscFunctionReturn(0); 190 } 191 192 PetscErrorCode MatDestroy_NormalHermitian(Mat N) { 193 Mat_Normal *Na = (Mat_Normal *)N->data; 194 195 PetscFunctionBegin; 196 PetscCall(MatDestroy(&Na->A)); 197 PetscCall(MatDestroy(&Na->D)); 198 PetscCall(VecDestroy(&Na->w)); 199 PetscCall(VecDestroy(&Na->left)); 200 PetscCall(VecDestroy(&Na->right)); 201 PetscCall(VecDestroy(&Na->leftwork)); 202 PetscCall(VecDestroy(&Na->rightwork)); 203 PetscCall(PetscFree(N->data)); 204 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMatHermitian_C", NULL)); 205 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL)); 206 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL)); 207 PetscFunctionReturn(0); 208 } 209 210 /* 211 Slow, nonscalable version 212 */ 213 PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v) { 214 Mat_Normal *Na = (Mat_Normal *)N->data; 215 Mat A = Na->A; 216 PetscInt i, j, rstart, rend, nnz; 217 const PetscInt *cols; 218 PetscScalar *diag, *work, *values; 219 const PetscScalar *mvalues; 220 221 PetscFunctionBegin; 222 PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work)); 223 PetscCall(PetscArrayzero(work, A->cmap->N)); 224 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 225 for (i = rstart; i < rend; i++) { 226 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues)); 227 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]); 228 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues)); 229 } 230 PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N))); 231 rstart = N->cmap->rstart; 232 rend = N->cmap->rend; 233 PetscCall(VecGetArray(v, &values)); 234 PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart)); 235 PetscCall(VecRestoreArray(v, &values)); 236 PetscCall(PetscFree2(diag, work)); 237 PetscCall(VecScale(v, Na->scale)); 238 PetscFunctionReturn(0); 239 } 240 241 PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D) { 242 Mat_Normal *Na = (Mat_Normal *)N->data; 243 Mat M, A = Na->A; 244 245 PetscFunctionBegin; 246 PetscCall(MatGetDiagonalBlock(A, &M)); 247 PetscCall(MatCreateNormalHermitian(M, &Na->D)); 248 *D = Na->D; 249 PetscFunctionReturn(0); 250 } 251 252 PetscErrorCode MatNormalGetMat_NormalHermitian(Mat A, Mat *M) { 253 Mat_Normal *Aa = (Mat_Normal *)A->data; 254 255 PetscFunctionBegin; 256 *M = Aa->A; 257 PetscFunctionReturn(0); 258 } 259 260 /*@ 261 MatNormalHermitianGetMat - Gets the Mat object stored inside a MATNORMALHERMITIAN 262 263 Logically collective on Mat 264 265 Input Parameter: 266 . A - the MATNORMALHERMITIAN matrix 267 268 Output Parameter: 269 . M - the matrix object stored inside A 270 271 Level: intermediate 272 273 .seealso: `MatCreateNormalHermitian()` 274 275 @*/ 276 PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M) { 277 PetscFunctionBegin; 278 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 279 PetscValidType(A, 1); 280 PetscValidPointer(M, 2); 281 PetscUseMethod(A, "MatNormalGetMatHermitian_C", (Mat, Mat *), (A, M)); 282 PetscFunctionReturn(0); 283 } 284 285 PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) { 286 Mat_Normal *Aa = (Mat_Normal *)A->data; 287 Mat B, conjugate; 288 PetscInt m, n, M, N; 289 290 PetscFunctionBegin; 291 PetscCall(MatGetSize(A, &M, &N)); 292 PetscCall(MatGetLocalSize(A, &m, &n)); 293 if (reuse == MAT_REUSE_MATRIX) { 294 B = *newmat; 295 PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B)); 296 } else { 297 PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B)); 298 PetscCall(MatProductSetType(B, MATPRODUCT_AtB)); 299 PetscCall(MatProductSetFromOptions(B)); 300 PetscCall(MatProductSymbolic(B)); 301 PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE)); 302 } 303 if (PetscDefined(USE_COMPLEX)) { 304 PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate)); 305 PetscCall(MatConjugate(conjugate)); 306 PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B)); 307 } 308 PetscCall(MatProductNumeric(B)); 309 if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate)); 310 if (reuse == MAT_INPLACE_MATRIX) { 311 PetscCall(MatHeaderReplace(A, &B)); 312 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 313 PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat)); 314 PetscFunctionReturn(0); 315 } 316 317 /*@ 318 MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A. 319 320 Collective on Mat 321 322 Input Parameter: 323 . A - the (possibly rectangular complex) matrix 324 325 Output Parameter: 326 . N - the matrix that represents (A*)'*A 327 328 Level: intermediate 329 330 Notes: 331 The product (A*)'*A is NOT actually formed! Rather the new matrix 332 object performs the matrix-vector product by first multiplying by 333 A and then (A*)' 334 @*/ 335 PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N) { 336 PetscInt m, n; 337 Mat_Normal *Na; 338 VecType vtype; 339 340 PetscFunctionBegin; 341 PetscCall(MatGetLocalSize(A, &m, &n)); 342 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 343 PetscCall(MatSetSizes(*N, n, n, PETSC_DECIDE, PETSC_DECIDE)); 344 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN)); 345 PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap)); 346 PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap)); 347 348 PetscCall(PetscNewLog(*N, &Na)); 349 (*N)->data = (void *)Na; 350 PetscCall(PetscObjectReference((PetscObject)A)); 351 Na->A = A; 352 Na->scale = 1.0; 353 354 PetscCall(MatCreateVecs(A, NULL, &Na->w)); 355 356 (*N)->ops->destroy = MatDestroy_NormalHermitian; 357 (*N)->ops->mult = MatMult_NormalHermitian; 358 (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal; 359 (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal; 360 (*N)->ops->multadd = MatMultHermitianAdd_Normal; 361 (*N)->ops->getdiagonal = MatGetDiagonal_NormalHermitian; 362 (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_NormalHermitian; 363 (*N)->ops->scale = MatScale_NormalHermitian; 364 (*N)->ops->diagonalscale = MatDiagonalScale_NormalHermitian; 365 (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian; 366 (*N)->ops->permute = MatPermute_NormalHermitian; 367 (*N)->ops->duplicate = MatDuplicate_NormalHermitian; 368 (*N)->ops->copy = MatCopy_NormalHermitian; 369 (*N)->assembled = PETSC_TRUE; 370 (*N)->preallocated = PETSC_TRUE; 371 372 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMatHermitian_C", MatNormalGetMat_NormalHermitian)); 373 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ)); 374 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ)); 375 PetscCall(MatSetOption(*N, MAT_HERMITIAN, PETSC_TRUE)); 376 PetscCall(MatGetVecType(A, &vtype)); 377 PetscCall(MatSetVecType(*N, vtype)); 378 #if defined(PETSC_HAVE_DEVICE) 379 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 380 #endif 381 PetscFunctionReturn(0); 382 } 383