1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; 6 Vec w,left,right,leftwork,rightwork; 7 PetscScalar scale; 8 } Mat_Normal; 9 10 PetscErrorCode MatScale_Normal(Mat inA,PetscScalar scale) 11 { 12 Mat_Normal *a = (Mat_Normal*)inA->data; 13 14 PetscFunctionBegin; 15 a->scale *= scale; 16 PetscFunctionReturn(0); 17 } 18 19 PetscErrorCode MatDiagonalScale_Normal(Mat inA,Vec left,Vec right) 20 { 21 Mat_Normal *a = (Mat_Normal*)inA->data; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 if (left) { 26 if (!a->left) { 27 ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr); 28 ierr = VecCopy(left,a->left);CHKERRQ(ierr); 29 } else { 30 ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr); 31 } 32 } 33 if (right) { 34 if (!a->right) { 35 ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr); 36 ierr = VecCopy(right,a->right);CHKERRQ(ierr); 37 } else { 38 ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr); 39 } 40 } 41 PetscFunctionReturn(0); 42 } 43 44 PetscErrorCode MatIncreaseOverlap_Normal(Mat A,PetscInt is_max,IS is[],PetscInt ov) 45 { 46 Mat_Normal *a = (Mat_Normal*)A->data; 47 Mat pattern; 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 PetscCheckFalse(ov < 0,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified"); 52 ierr = MatProductCreate(a->A,a->A,NULL,&pattern);CHKERRQ(ierr); 53 ierr = MatProductSetType(pattern,MATPRODUCT_AtB);CHKERRQ(ierr); 54 ierr = MatProductSetFromOptions(pattern);CHKERRQ(ierr); 55 ierr = MatProductSymbolic(pattern);CHKERRQ(ierr); 56 ierr = MatIncreaseOverlap(pattern,is_max,is,ov);CHKERRQ(ierr); 57 ierr = MatDestroy(&pattern);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 PetscErrorCode MatCreateSubMatrices_Normal(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[]) 62 { 63 Mat_Normal *a = (Mat_Normal*)mat->data; 64 Mat B = a->A, *suba; 65 IS *row; 66 PetscInt M; 67 PetscErrorCode ierr; 68 69 PetscFunctionBegin; 70 PetscCheckFalse(a->left || a->right || irow != icol,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Not implemented"); 71 if (scall != MAT_REUSE_MATRIX) { 72 ierr = PetscCalloc1(n,submat);CHKERRQ(ierr); 73 } 74 ierr = MatGetSize(B,&M,NULL);CHKERRQ(ierr); 75 ierr = PetscMalloc1(n,&row);CHKERRQ(ierr); 76 ierr = ISCreateStride(PETSC_COMM_SELF,M,0,1,&row[0]);CHKERRQ(ierr); 77 ierr = ISSetIdentity(row[0]);CHKERRQ(ierr); 78 for (M = 1; M < n; ++M) row[M] = row[0]; 79 ierr = MatCreateSubMatrices(B,n,row,icol,MAT_INITIAL_MATRIX,&suba);CHKERRQ(ierr); 80 for (M = 0; M < n; ++M) { 81 ierr = MatCreateNormal(suba[M],*submat+M);CHKERRQ(ierr); 82 ((Mat_Normal*)(*submat)[M]->data)->scale = a->scale; 83 } 84 ierr = ISDestroy(&row[0]);CHKERRQ(ierr); 85 ierr = PetscFree(row);CHKERRQ(ierr); 86 ierr = MatDestroySubMatrices(n,&suba);CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 PetscErrorCode MatPermute_Normal(Mat A,IS rowp,IS colp,Mat *B) 91 { 92 Mat_Normal *a = (Mat_Normal*)A->data; 93 Mat C,Aa = a->A; 94 IS row; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 PetscCheckFalse(rowp != colp,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Row permutation and column permutation must be the same"); 99 ierr = ISCreateStride(PetscObjectComm((PetscObject)Aa),Aa->rmap->n,Aa->rmap->rstart,1,&row);CHKERRQ(ierr); 100 ierr = ISSetIdentity(row);CHKERRQ(ierr); 101 ierr = MatPermute(Aa,row,colp,&C);CHKERRQ(ierr); 102 ierr = ISDestroy(&row);CHKERRQ(ierr); 103 ierr = MatCreateNormal(C,B);CHKERRQ(ierr); 104 ierr = MatDestroy(&C);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 108 PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B) 109 { 110 Mat_Normal *a = (Mat_Normal*)A->data; 111 Mat C; 112 PetscErrorCode ierr; 113 114 PetscFunctionBegin; 115 PetscCheckFalse(a->left || a->right,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented"); 116 ierr = MatDuplicate(a->A,op,&C);CHKERRQ(ierr); 117 ierr = MatCreateNormal(C,B);CHKERRQ(ierr); 118 ierr = MatDestroy(&C);CHKERRQ(ierr); 119 if (op == MAT_COPY_VALUES) ((Mat_Normal*)(*B)->data)->scale = a->scale; 120 PetscFunctionReturn(0); 121 } 122 123 PetscErrorCode MatCopy_Normal(Mat A,Mat B,MatStructure str) 124 { 125 Mat_Normal *a = (Mat_Normal*)A->data,*b = (Mat_Normal*)B->data; 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 PetscCheckFalse(a->left || a->right,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented"); 130 ierr = MatCopy(a->A,b->A,str);CHKERRQ(ierr); 131 b->scale = a->scale; 132 ierr = VecDestroy(&b->left);CHKERRQ(ierr); 133 ierr = VecDestroy(&b->right);CHKERRQ(ierr); 134 ierr = VecDestroy(&b->leftwork);CHKERRQ(ierr); 135 ierr = VecDestroy(&b->rightwork);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y) 140 { 141 Mat_Normal *Na = (Mat_Normal*)N->data; 142 PetscErrorCode ierr; 143 Vec in; 144 145 PetscFunctionBegin; 146 in = x; 147 if (Na->right) { 148 if (!Na->rightwork) { 149 ierr = VecDuplicate(Na->right,&Na->rightwork);CHKERRQ(ierr); 150 } 151 ierr = VecPointwiseMult(Na->rightwork,Na->right,in);CHKERRQ(ierr); 152 in = Na->rightwork; 153 } 154 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr); 155 ierr = MatMultTranspose(Na->A,Na->w,y);CHKERRQ(ierr); 156 if (Na->left) { 157 ierr = VecPointwiseMult(y,Na->left,y);CHKERRQ(ierr); 158 } 159 ierr = VecScale(y,Na->scale);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3) 164 { 165 Mat_Normal *Na = (Mat_Normal*)N->data; 166 PetscErrorCode ierr; 167 Vec in; 168 169 PetscFunctionBegin; 170 in = v1; 171 if (Na->right) { 172 if (!Na->rightwork) { 173 ierr = VecDuplicate(Na->right,&Na->rightwork);CHKERRQ(ierr); 174 } 175 ierr = VecPointwiseMult(Na->rightwork,Na->right,in);CHKERRQ(ierr); 176 in = Na->rightwork; 177 } 178 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr); 179 ierr = VecScale(Na->w,Na->scale);CHKERRQ(ierr); 180 if (Na->left) { 181 ierr = MatMultTranspose(Na->A,Na->w,v3);CHKERRQ(ierr); 182 ierr = VecPointwiseMult(v3,Na->left,v3);CHKERRQ(ierr); 183 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 184 } else { 185 ierr = MatMultTransposeAdd(Na->A,Na->w,v2,v3);CHKERRQ(ierr); 186 } 187 PetscFunctionReturn(0); 188 } 189 190 PetscErrorCode MatMultTranspose_Normal(Mat N,Vec x,Vec y) 191 { 192 Mat_Normal *Na = (Mat_Normal*)N->data; 193 PetscErrorCode ierr; 194 Vec in; 195 196 PetscFunctionBegin; 197 in = x; 198 if (Na->left) { 199 if (!Na->leftwork) { 200 ierr = VecDuplicate(Na->left,&Na->leftwork);CHKERRQ(ierr); 201 } 202 ierr = VecPointwiseMult(Na->leftwork,Na->left,in);CHKERRQ(ierr); 203 in = Na->leftwork; 204 } 205 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr); 206 ierr = MatMultTranspose(Na->A,Na->w,y);CHKERRQ(ierr); 207 if (Na->right) { 208 ierr = VecPointwiseMult(y,Na->right,y);CHKERRQ(ierr); 209 } 210 ierr = VecScale(y,Na->scale);CHKERRQ(ierr); 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode MatMultTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3) 215 { 216 Mat_Normal *Na = (Mat_Normal*)N->data; 217 PetscErrorCode ierr; 218 Vec in; 219 220 PetscFunctionBegin; 221 in = v1; 222 if (Na->left) { 223 if (!Na->leftwork) { 224 ierr = VecDuplicate(Na->left,&Na->leftwork);CHKERRQ(ierr); 225 } 226 ierr = VecPointwiseMult(Na->leftwork,Na->left,in);CHKERRQ(ierr); 227 in = Na->leftwork; 228 } 229 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr); 230 ierr = VecScale(Na->w,Na->scale);CHKERRQ(ierr); 231 if (Na->right) { 232 ierr = MatMultTranspose(Na->A,Na->w,v3);CHKERRQ(ierr); 233 ierr = VecPointwiseMult(v3,Na->right,v3);CHKERRQ(ierr); 234 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 235 } else { 236 ierr = MatMultTransposeAdd(Na->A,Na->w,v2,v3);CHKERRQ(ierr); 237 } 238 PetscFunctionReturn(0); 239 } 240 241 PetscErrorCode MatDestroy_Normal(Mat N) 242 { 243 Mat_Normal *Na = (Mat_Normal*)N->data; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 ierr = MatDestroy(&Na->A);CHKERRQ(ierr); 248 ierr = VecDestroy(&Na->w);CHKERRQ(ierr); 249 ierr = VecDestroy(&Na->left);CHKERRQ(ierr); 250 ierr = VecDestroy(&Na->right);CHKERRQ(ierr); 251 ierr = VecDestroy(&Na->leftwork);CHKERRQ(ierr); 252 ierr = VecDestroy(&Na->rightwork);CHKERRQ(ierr); 253 ierr = PetscFree(N->data);CHKERRQ(ierr); 254 ierr = PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMat_C",NULL);CHKERRQ(ierr); 255 ierr = PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_seqaij_C",NULL);CHKERRQ(ierr); 256 ierr = PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_mpiaij_C",NULL);CHKERRQ(ierr); 257 ierr = PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_seqdense_C",NULL);CHKERRQ(ierr); 258 ierr = PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_mpidense_C",NULL);CHKERRQ(ierr); 259 ierr = PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_dense_C",NULL);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 263 /* 264 Slow, nonscalable version 265 */ 266 PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v) 267 { 268 Mat_Normal *Na = (Mat_Normal*)N->data; 269 Mat A = Na->A; 270 PetscErrorCode ierr; 271 PetscInt i,j,rstart,rend,nnz; 272 const PetscInt *cols; 273 PetscScalar *diag,*work,*values; 274 const PetscScalar *mvalues; 275 276 PetscFunctionBegin; 277 ierr = PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);CHKERRQ(ierr); 278 ierr = PetscArrayzero(work,A->cmap->N);CHKERRQ(ierr); 279 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 280 for (i=rstart; i<rend; i++) { 281 ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 282 for (j=0; j<nnz; j++) { 283 work[cols[j]] += mvalues[j]*mvalues[j]; 284 } 285 ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 286 } 287 ierr = MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRMPI(ierr); 288 rstart = N->cmap->rstart; 289 rend = N->cmap->rend; 290 ierr = VecGetArray(v,&values);CHKERRQ(ierr); 291 ierr = PetscArraycpy(values,diag+rstart,rend-rstart);CHKERRQ(ierr); 292 ierr = VecRestoreArray(v,&values);CHKERRQ(ierr); 293 ierr = PetscFree2(diag,work);CHKERRQ(ierr); 294 ierr = VecScale(v,Na->scale);CHKERRQ(ierr); 295 PetscFunctionReturn(0); 296 } 297 298 PetscErrorCode MatNormalGetMat_Normal(Mat A,Mat *M) 299 { 300 Mat_Normal *Aa = (Mat_Normal*)A->data; 301 302 PetscFunctionBegin; 303 *M = Aa->A; 304 PetscFunctionReturn(0); 305 } 306 307 /*@ 308 MatNormalGetMat - Gets the Mat object stored inside a MATNORMAL 309 310 Logically collective on Mat 311 312 Input Parameter: 313 . A - the MATNORMAL matrix 314 315 Output Parameter: 316 . M - the matrix object stored inside A 317 318 Level: intermediate 319 320 .seealso: MatCreateNormal() 321 322 @*/ 323 PetscErrorCode MatNormalGetMat(Mat A,Mat *M) 324 { 325 PetscErrorCode ierr; 326 327 PetscFunctionBegin; 328 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 329 PetscValidType(A,1); 330 PetscValidPointer(M,2); 331 ierr = PetscUseMethod(A,"MatNormalGetMat_C",(Mat,Mat*),(A,M));CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 PetscErrorCode MatConvert_Normal_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 336 { 337 Mat_Normal *Aa = (Mat_Normal*)A->data; 338 Mat B; 339 PetscInt m,n,M,N; 340 PetscErrorCode ierr; 341 342 PetscFunctionBegin; 343 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 344 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 345 if (reuse == MAT_REUSE_MATRIX) { 346 B = *newmat; 347 ierr = MatProductReplaceMats(Aa->A,Aa->A,NULL,B);CHKERRQ(ierr); 348 } else { 349 ierr = MatProductCreate(Aa->A,Aa->A,NULL,&B);CHKERRQ(ierr); 350 ierr = MatProductSetType(B,MATPRODUCT_AtB);CHKERRQ(ierr); 351 ierr = MatProductSetFromOptions(B);CHKERRQ(ierr); 352 ierr = MatProductSymbolic(B);CHKERRQ(ierr); 353 ierr = MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 354 } 355 ierr = MatProductNumeric(B);CHKERRQ(ierr); 356 if (reuse == MAT_INPLACE_MATRIX) { 357 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 358 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 359 ierr = MatConvert(*newmat,MATAIJ,MAT_INPLACE_MATRIX,newmat);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 typedef struct { 364 Mat work[2]; 365 } Normal_Dense; 366 367 PetscErrorCode MatProductNumeric_Normal_Dense(Mat C) 368 { 369 Mat A,B; 370 Normal_Dense *contents; 371 Mat_Normal *a; 372 PetscScalar *array; 373 PetscErrorCode ierr; 374 375 PetscFunctionBegin; 376 MatCheckProduct(C,3); 377 A = C->product->A; 378 a = (Mat_Normal*)A->data; 379 B = C->product->B; 380 contents = (Normal_Dense*)C->product->data; 381 PetscCheckFalse(!contents,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty"); 382 if (a->right) { 383 ierr = MatCopy(B,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 384 ierr = MatDiagonalScale(C,a->right,NULL);CHKERRQ(ierr); 385 } 386 ierr = MatProductNumeric(contents->work[0]);CHKERRQ(ierr); 387 ierr = MatDenseGetArrayWrite(C,&array);CHKERRQ(ierr); 388 ierr = MatDensePlaceArray(contents->work[1],array);CHKERRQ(ierr); 389 ierr = MatProductNumeric(contents->work[1]);CHKERRQ(ierr); 390 ierr = MatDenseRestoreArrayWrite(C,&array);CHKERRQ(ierr); 391 ierr = MatDenseResetArray(contents->work[1]);CHKERRQ(ierr); 392 ierr = MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 393 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 394 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 395 ierr = MatScale(C,a->scale);CHKERRQ(ierr); 396 PetscFunctionReturn(0); 397 } 398 399 PetscErrorCode MatNormal_DenseDestroy(void *ctx) 400 { 401 Normal_Dense *contents = (Normal_Dense*)ctx; 402 PetscErrorCode ierr; 403 404 PetscFunctionBegin; 405 ierr = MatDestroy(contents->work);CHKERRQ(ierr); 406 ierr = MatDestroy(contents->work+1);CHKERRQ(ierr); 407 ierr = PetscFree(contents);CHKERRQ(ierr); 408 PetscFunctionReturn(0); 409 } 410 411 PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C) 412 { 413 Mat A,B; 414 Normal_Dense *contents = NULL; 415 Mat_Normal *a; 416 PetscScalar *array; 417 PetscInt n,N,m,M; 418 PetscErrorCode ierr; 419 420 PetscFunctionBegin; 421 MatCheckProduct(C,4); 422 PetscCheckFalse(C->product->data,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty"); 423 A = C->product->A; 424 a = (Mat_Normal*)A->data; 425 PetscCheckFalse(a->left,PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Not implemented"); 426 B = C->product->B; 427 ierr = MatGetLocalSize(C,&m,&n);CHKERRQ(ierr); 428 ierr = MatGetSize(C,&M,&N);CHKERRQ(ierr); 429 if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) { 430 ierr = MatGetLocalSize(B,NULL,&n);CHKERRQ(ierr); 431 ierr = MatGetSize(B,NULL,&N);CHKERRQ(ierr); 432 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 433 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 434 ierr = MatSetSizes(C,m,n,M,N);CHKERRQ(ierr); 435 } 436 ierr = MatSetType(C,((PetscObject)B)->type_name);CHKERRQ(ierr); 437 ierr = MatSetUp(C);CHKERRQ(ierr); 438 ierr = PetscNew(&contents);CHKERRQ(ierr); 439 C->product->data = contents; 440 C->product->destroy = MatNormal_DenseDestroy; 441 if (a->right) { 442 ierr = MatProductCreate(a->A,C,NULL,contents->work);CHKERRQ(ierr); 443 } else { 444 ierr = MatProductCreate(a->A,B,NULL,contents->work);CHKERRQ(ierr); 445 } 446 ierr = MatProductSetType(contents->work[0],MATPRODUCT_AB);CHKERRQ(ierr); 447 ierr = MatProductSetFromOptions(contents->work[0]);CHKERRQ(ierr); 448 ierr = MatProductSymbolic(contents->work[0]);CHKERRQ(ierr); 449 ierr = MatProductCreate(a->A,contents->work[0],NULL,contents->work+1);CHKERRQ(ierr); 450 ierr = MatProductSetType(contents->work[1],MATPRODUCT_AtB);CHKERRQ(ierr); 451 ierr = MatProductSetFromOptions(contents->work[1]);CHKERRQ(ierr); 452 ierr = MatProductSymbolic(contents->work[1]);CHKERRQ(ierr); 453 ierr = MatDenseGetArrayWrite(C,&array);CHKERRQ(ierr); 454 ierr = MatSeqDenseSetPreallocation(contents->work[1],array);CHKERRQ(ierr); 455 ierr = MatMPIDenseSetPreallocation(contents->work[1],array);CHKERRQ(ierr); 456 ierr = MatDenseRestoreArrayWrite(C,&array);CHKERRQ(ierr); 457 C->ops->productnumeric = MatProductNumeric_Normal_Dense; 458 PetscFunctionReturn(0); 459 } 460 461 PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C) 462 { 463 PetscFunctionBegin; 464 C->ops->productsymbolic = MatProductSymbolic_Normal_Dense; 465 PetscFunctionReturn(0); 466 } 467 468 PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C) 469 { 470 Mat_Product *product = C->product; 471 PetscErrorCode ierr; 472 473 PetscFunctionBegin; 474 if (product->type == MATPRODUCT_AB) { 475 ierr = MatProductSetFromOptions_Normal_Dense_AB(C);CHKERRQ(ierr); 476 } 477 PetscFunctionReturn(0); 478 } 479 480 /*@ 481 MatCreateNormal - Creates a new matrix object that behaves like A'*A. 482 483 Collective on Mat 484 485 Input Parameter: 486 . A - the (possibly rectangular) matrix 487 488 Output Parameter: 489 . N - the matrix that represents A'*A 490 491 Level: intermediate 492 493 Notes: 494 The product A'*A is NOT actually formed! Rather the new matrix 495 object performs the matrix-vector product by first multiplying by 496 A and then A' 497 @*/ 498 PetscErrorCode MatCreateNormal(Mat A,Mat *N) 499 { 500 PetscErrorCode ierr; 501 PetscInt n,nn; 502 Mat_Normal *Na; 503 VecType vtype; 504 505 PetscFunctionBegin; 506 ierr = MatGetSize(A,NULL,&nn);CHKERRQ(ierr); 507 ierr = MatGetLocalSize(A,NULL,&n);CHKERRQ(ierr); 508 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 509 ierr = MatSetSizes(*N,n,n,nn,nn);CHKERRQ(ierr); 510 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);CHKERRQ(ierr); 511 ierr = PetscLayoutReference(A->cmap,&(*N)->rmap);CHKERRQ(ierr); 512 ierr = PetscLayoutReference(A->cmap,&(*N)->cmap);CHKERRQ(ierr); 513 514 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 515 (*N)->data = (void*) Na; 516 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 517 Na->A = A; 518 Na->scale = 1.0; 519 520 ierr = MatCreateVecs(A,NULL,&Na->w);CHKERRQ(ierr); 521 522 (*N)->ops->destroy = MatDestroy_Normal; 523 (*N)->ops->mult = MatMult_Normal; 524 (*N)->ops->multtranspose = MatMultTranspose_Normal; 525 (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal; 526 (*N)->ops->multadd = MatMultAdd_Normal; 527 (*N)->ops->getdiagonal = MatGetDiagonal_Normal; 528 (*N)->ops->scale = MatScale_Normal; 529 (*N)->ops->diagonalscale = MatDiagonalScale_Normal; 530 (*N)->ops->increaseoverlap = MatIncreaseOverlap_Normal; 531 (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal; 532 (*N)->ops->permute = MatPermute_Normal; 533 (*N)->ops->duplicate = MatDuplicate_Normal; 534 (*N)->ops->copy = MatCopy_Normal; 535 (*N)->assembled = PETSC_TRUE; 536 (*N)->preallocated = PETSC_TRUE; 537 538 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMat_C",MatNormalGetMat_Normal);CHKERRQ(ierr); 539 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_seqaij_C",MatConvert_Normal_AIJ);CHKERRQ(ierr); 540 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_mpiaij_C",MatConvert_Normal_AIJ);CHKERRQ(ierr); 541 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_seqdense_C",MatProductSetFromOptions_Normal_Dense);CHKERRQ(ierr); 542 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_mpidense_C",MatProductSetFromOptions_Normal_Dense);CHKERRQ(ierr); 543 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_dense_C",MatProductSetFromOptions_Normal_Dense);CHKERRQ(ierr); 544 ierr = MatSetOption(*N,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 545 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr); 546 ierr = MatSetVecType(*N,vtype);CHKERRQ(ierr); 547 #if defined(PETSC_HAVE_DEVICE) 548 ierr = MatBindToCPU(*N,A->boundtocpu);CHKERRQ(ierr); 549 #endif 550 PetscFunctionReturn(0); 551 } 552