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 MatMult_NormalHermitian(Mat N,Vec x,Vec y) 45 { 46 Mat_Normal *Na = (Mat_Normal*)N->data; 47 Vec in; 48 49 PetscFunctionBegin; 50 in = x; 51 if (Na->right) { 52 if (!Na->rightwork) { 53 PetscCall(VecDuplicate(Na->right,&Na->rightwork)); 54 } 55 PetscCall(VecPointwiseMult(Na->rightwork,Na->right,in)); 56 in = Na->rightwork; 57 } 58 PetscCall(MatMult(Na->A,in,Na->w)); 59 PetscCall(MatMultHermitianTranspose(Na->A,Na->w,y)); 60 if (Na->left) { 61 PetscCall(VecPointwiseMult(y,Na->left,y)); 62 } 63 PetscCall(VecScale(y,Na->scale)); 64 PetscFunctionReturn(0); 65 } 66 67 PetscErrorCode MatMultHermitianAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3) 68 { 69 Mat_Normal *Na = (Mat_Normal*)N->data; 70 Vec in; 71 72 PetscFunctionBegin; 73 in = v1; 74 if (Na->right) { 75 if (!Na->rightwork) { 76 PetscCall(VecDuplicate(Na->right,&Na->rightwork)); 77 } 78 PetscCall(VecPointwiseMult(Na->rightwork,Na->right,in)); 79 in = Na->rightwork; 80 } 81 PetscCall(MatMult(Na->A,in,Na->w)); 82 PetscCall(VecScale(Na->w,Na->scale)); 83 if (Na->left) { 84 PetscCall(MatMultHermitianTranspose(Na->A,Na->w,v3)); 85 PetscCall(VecPointwiseMult(v3,Na->left,v3)); 86 PetscCall(VecAXPY(v3,1.0,v2)); 87 } else { 88 PetscCall(MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3)); 89 } 90 PetscFunctionReturn(0); 91 } 92 93 PetscErrorCode MatMultHermitianTranspose_Normal(Mat N,Vec x,Vec y) 94 { 95 Mat_Normal *Na = (Mat_Normal*)N->data; 96 Vec in; 97 98 PetscFunctionBegin; 99 in = x; 100 if (Na->left) { 101 if (!Na->leftwork) { 102 PetscCall(VecDuplicate(Na->left,&Na->leftwork)); 103 } 104 PetscCall(VecPointwiseMult(Na->leftwork,Na->left,in)); 105 in = Na->leftwork; 106 } 107 PetscCall(MatMult(Na->A,in,Na->w)); 108 PetscCall(MatMultHermitianTranspose(Na->A,Na->w,y)); 109 if (Na->right) { 110 PetscCall(VecPointwiseMult(y,Na->right,y)); 111 } 112 PetscCall(VecScale(y,Na->scale)); 113 PetscFunctionReturn(0); 114 } 115 116 PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3) 117 { 118 Mat_Normal *Na = (Mat_Normal*)N->data; 119 Vec in; 120 121 PetscFunctionBegin; 122 in = v1; 123 if (Na->left) { 124 if (!Na->leftwork) { 125 PetscCall(VecDuplicate(Na->left,&Na->leftwork)); 126 } 127 PetscCall(VecPointwiseMult(Na->leftwork,Na->left,in)); 128 in = Na->leftwork; 129 } 130 PetscCall(MatMult(Na->A,in,Na->w)); 131 PetscCall(VecScale(Na->w,Na->scale)); 132 if (Na->right) { 133 PetscCall(MatMultHermitianTranspose(Na->A,Na->w,v3)); 134 PetscCall(VecPointwiseMult(v3,Na->right,v3)); 135 PetscCall(VecAXPY(v3,1.0,v2)); 136 } else { 137 PetscCall(MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3)); 138 } 139 PetscFunctionReturn(0); 140 } 141 142 PetscErrorCode MatDestroy_NormalHermitian(Mat N) 143 { 144 Mat_Normal *Na = (Mat_Normal*)N->data; 145 146 PetscFunctionBegin; 147 PetscCall(MatDestroy(&Na->A)); 148 PetscCall(MatDestroy(&Na->D)); 149 PetscCall(VecDestroy(&Na->w)); 150 PetscCall(VecDestroy(&Na->left)); 151 PetscCall(VecDestroy(&Na->right)); 152 PetscCall(VecDestroy(&Na->leftwork)); 153 PetscCall(VecDestroy(&Na->rightwork)); 154 PetscCall(PetscFree(N->data)); 155 PetscCall(PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMatHermitian_C",NULL)); 156 PetscFunctionReturn(0); 157 } 158 159 /* 160 Slow, nonscalable version 161 */ 162 PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N,Vec v) 163 { 164 Mat_Normal *Na = (Mat_Normal*)N->data; 165 Mat A = Na->A; 166 PetscInt i,j,rstart,rend,nnz; 167 const PetscInt *cols; 168 PetscScalar *diag,*work,*values; 169 const PetscScalar *mvalues; 170 171 PetscFunctionBegin; 172 PetscCall(PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work)); 173 PetscCall(PetscArrayzero(work,A->cmap->N)); 174 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 175 for (i=rstart; i<rend; i++) { 176 PetscCall(MatGetRow(A,i,&nnz,&cols,&mvalues)); 177 for (j=0; j<nnz; j++) { 178 work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]); 179 } 180 PetscCall(MatRestoreRow(A,i,&nnz,&cols,&mvalues)); 181 } 182 PetscCall(MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N))); 183 rstart = N->cmap->rstart; 184 rend = N->cmap->rend; 185 PetscCall(VecGetArray(v,&values)); 186 PetscCall(PetscArraycpy(values,diag+rstart,rend-rstart)); 187 PetscCall(VecRestoreArray(v,&values)); 188 PetscCall(PetscFree2(diag,work)); 189 PetscCall(VecScale(v,Na->scale)); 190 PetscFunctionReturn(0); 191 } 192 193 PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N,Mat *D) 194 { 195 Mat_Normal *Na = (Mat_Normal*)N->data; 196 Mat M,A = Na->A; 197 198 PetscFunctionBegin; 199 PetscCall(MatGetDiagonalBlock(A,&M)); 200 PetscCall(MatCreateNormalHermitian(M,&Na->D)); 201 *D = Na->D; 202 PetscFunctionReturn(0); 203 } 204 205 PetscErrorCode MatNormalGetMat_NormalHermitian(Mat A,Mat *M) 206 { 207 Mat_Normal *Aa = (Mat_Normal*)A->data; 208 209 PetscFunctionBegin; 210 *M = Aa->A; 211 PetscFunctionReturn(0); 212 } 213 214 /*@ 215 MatNormalHermitianGetMat - Gets the Mat object stored inside a MATNORMALHERMITIAN 216 217 Logically collective on Mat 218 219 Input Parameter: 220 . A - the MATNORMALHERMITIAN matrix 221 222 Output Parameter: 223 . M - the matrix object stored inside A 224 225 Level: intermediate 226 227 .seealso: `MatCreateNormalHermitian()` 228 229 @*/ 230 PetscErrorCode MatNormalHermitianGetMat(Mat A,Mat *M) 231 { 232 PetscFunctionBegin; 233 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 234 PetscValidType(A,1); 235 PetscValidPointer(M,2); 236 PetscUseMethod(A,"MatNormalGetMatHermitian_C",(Mat,Mat*),(A,M)); 237 PetscFunctionReturn(0); 238 } 239 240 /*@ 241 MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A. 242 243 Collective on Mat 244 245 Input Parameter: 246 . A - the (possibly rectangular complex) matrix 247 248 Output Parameter: 249 . N - the matrix that represents (A*)'*A 250 251 Level: intermediate 252 253 Notes: 254 The product (A*)'*A is NOT actually formed! Rather the new matrix 255 object performs the matrix-vector product by first multiplying by 256 A and then (A*)' 257 @*/ 258 PetscErrorCode MatCreateNormalHermitian(Mat A,Mat *N) 259 { 260 PetscInt m,n; 261 Mat_Normal *Na; 262 VecType vtype; 263 264 PetscFunctionBegin; 265 PetscCall(MatGetLocalSize(A,&m,&n)); 266 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),N)); 267 PetscCall(MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE)); 268 PetscCall(PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN)); 269 PetscCall(PetscLayoutReference(A->cmap,&(*N)->rmap)); 270 PetscCall(PetscLayoutReference(A->cmap,&(*N)->cmap)); 271 272 PetscCall(PetscNewLog(*N,&Na)); 273 (*N)->data = (void*) Na; 274 PetscCall(PetscObjectReference((PetscObject)A)); 275 Na->A = A; 276 Na->scale = 1.0; 277 278 PetscCall(MatCreateVecs(A,NULL,&Na->w)); 279 280 (*N)->ops->destroy = MatDestroy_NormalHermitian; 281 (*N)->ops->mult = MatMult_NormalHermitian; 282 (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal; 283 (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal; 284 (*N)->ops->multadd = MatMultHermitianAdd_Normal; 285 (*N)->ops->getdiagonal = MatGetDiagonal_NormalHermitian; 286 (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_NormalHermitian; 287 (*N)->ops->scale = MatScale_NormalHermitian; 288 (*N)->ops->diagonalscale = MatDiagonalScale_NormalHermitian; 289 (*N)->assembled = PETSC_TRUE; 290 (*N)->preallocated = PETSC_TRUE; 291 292 PetscCall(PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMatHermitian_C",MatNormalGetMat_NormalHermitian)); 293 PetscCall(MatSetOption(*N,MAT_HERMITIAN,PETSC_TRUE)); 294 PetscCall(MatGetVecType(A,&vtype)); 295 PetscCall(MatSetVecType(*N,vtype)); 296 #if defined(PETSC_HAVE_DEVICE) 297 PetscCall(MatBindToCPU(*N,A->boundtocpu)); 298 #endif 299 PetscFunctionReturn(0); 300 } 301