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