1 #include <petscmat.h> /*I "mat.h" I*/ 2 3 typedef struct{ 4 Mat A; 5 Vec D1; 6 Vec D2; 7 Vec W; 8 Vec W2; 9 Vec ADADiag; 10 PetscInt GotDiag; 11 } _p_TaoMatADACtx; 12 typedef _p_TaoMatADACtx* TaoMatADACtx; 13 14 extern PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*); 15 16 static PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y) 17 { 18 TaoMatADACtx ctx; 19 PetscReal one = 1.0; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 24 ierr = MatMult(ctx->A,a,ctx->W);CHKERRQ(ierr); 25 if (ctx->D1){ 26 ierr = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(ierr); 27 } 28 ierr = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(ierr); 29 if (ctx->D2){ 30 ierr = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(ierr); 31 ierr = VecAXPY(y, one, ctx->W2);CHKERRQ(ierr); 32 } 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y) 37 { 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = MatMult_ADA(mat,a,y);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 static PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M) 46 { 47 TaoMatADACtx ctx; 48 PetscReal zero=0.0,one = 1.0; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr); 53 if (!ctx->D2){ 54 ierr = VecDuplicate(D,&ctx->D2);CHKERRQ(ierr); 55 ierr = VecSet(ctx->D2, zero);CHKERRQ(ierr); 56 } 57 ierr = VecAXPY(ctx->D2, one, D);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 static PetscErrorCode MatDestroy_ADA(Mat mat) 62 { 63 PetscErrorCode ierr; 64 TaoMatADACtx ctx; 65 66 PetscFunctionBegin; 67 ierr=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 68 ierr=VecDestroy(&ctx->W);CHKERRQ(ierr); 69 ierr=VecDestroy(&ctx->W2);CHKERRQ(ierr); 70 ierr=VecDestroy(&ctx->ADADiag);CHKERRQ(ierr); 71 ierr=MatDestroy(&ctx->A);CHKERRQ(ierr); 72 ierr=VecDestroy(&ctx->D1);CHKERRQ(ierr); 73 ierr=VecDestroy(&ctx->D2);CHKERRQ(ierr); 74 ierr = PetscFree(ctx);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer) 79 { 80 PetscFunctionBegin; 81 PetscFunctionReturn(0); 82 } 83 84 static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a) 85 { 86 PetscErrorCode ierr; 87 TaoMatADACtx ctx; 88 89 PetscFunctionBegin; 90 ierr = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(ierr); 91 ierr = VecShift(ctx->D2,a);CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M) 96 { 97 PetscErrorCode ierr; 98 TaoMatADACtx ctx; 99 Mat A2; 100 Vec D1b=NULL,D2b; 101 102 PetscFunctionBegin; 103 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 104 ierr = MatDuplicate(ctx->A,op,&A2);CHKERRQ(ierr); 105 if (ctx->D1){ 106 ierr = VecDuplicate(ctx->D1,&D1b);CHKERRQ(ierr); 107 ierr = VecCopy(ctx->D1,D1b);CHKERRQ(ierr); 108 } 109 ierr = VecDuplicate(ctx->D2,&D2b);CHKERRQ(ierr); 110 ierr = VecCopy(ctx->D2,D2b);CHKERRQ(ierr); 111 ierr = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(ierr); 112 if (ctx->D1){ 113 ierr = PetscObjectDereference((PetscObject)D1b);CHKERRQ(ierr); 114 } 115 ierr = PetscObjectDereference((PetscObject)D2b);CHKERRQ(ierr); 116 ierr = PetscObjectDereference((PetscObject)A2);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg) 121 { 122 PetscErrorCode ierr; 123 TaoMatADACtx ctx1,ctx2; 124 125 PetscFunctionBegin; 126 ierr = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(ierr); 127 ierr = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(ierr); 128 ierr = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(ierr); 129 if (*flg==PETSC_TRUE){ 130 ierr = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(ierr); 131 } 132 if (*flg==PETSC_TRUE){ 133 ierr = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(ierr); 134 } 135 PetscFunctionReturn(0); 136 } 137 138 static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a) 139 { 140 PetscErrorCode ierr; 141 TaoMatADACtx ctx; 142 143 PetscFunctionBegin; 144 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 145 ierr = VecScale(ctx->D1,a);CHKERRQ(ierr); 146 if (ctx->D2){ 147 ierr = VecScale(ctx->D2,a);CHKERRQ(ierr); 148 } 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode MatTranspose_ADA(Mat mat,MatReuse reuse,Mat *B) 153 { 154 PetscErrorCode ierr; 155 TaoMatADACtx ctx; 156 157 PetscFunctionBegin; 158 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 159 if (reuse == MAT_INITIAL_MATRIX){ 160 ierr = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(ierr); 161 } else if (reuse == MAT_REUSE_MATRIX){ 162 ierr = MatCopy(mat,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 163 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Does not support inplace transpose"); 164 PetscFunctionReturn(0); 165 } 166 167 PetscErrorCode MatADAComputeDiagonal(Mat mat) 168 { 169 PetscErrorCode ierr; 170 PetscInt i,m,n,low,high; 171 PetscScalar *dtemp,*dptr; 172 TaoMatADACtx ctx; 173 174 PetscFunctionBegin; 175 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 176 ierr = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(ierr); 177 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 178 179 ierr = PetscMalloc1(n,&dtemp );CHKERRQ(ierr); 180 181 for (i=0; i<n; i++){ 182 ierr = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(ierr); 183 ierr = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(ierr); 184 ierr = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 185 } 186 for (i=0; i<n; i++){ 187 ierr = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 188 } 189 190 ierr = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 191 for (i=low; i<high; i++){ 192 dptr[i-low]= dtemp[i]; 193 } 194 ierr = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 195 ierr = PetscFree(dtemp);CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 199 static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v) 200 { 201 PetscErrorCode ierr; 202 PetscReal one=1.0; 203 TaoMatADACtx ctx; 204 205 PetscFunctionBegin; 206 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 207 ierr = MatADAComputeDiagonal(mat);CHKERRQ(ierr); 208 ierr=VecCopy(ctx->ADADiag,v);CHKERRQ(ierr); 209 if (ctx->D2){ 210 ierr=VecAXPY(v, one, ctx->D2);CHKERRQ(ierr); 211 } 212 PetscFunctionReturn(0); 213 } 214 215 static PetscErrorCode MatCreateSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat) 216 { 217 PetscErrorCode ierr; 218 PetscInt low,high; 219 IS ISrow; 220 Vec D1,D2; 221 Mat Atemp; 222 TaoMatADACtx ctx; 223 PetscBool isequal; 224 225 PetscFunctionBegin; 226 ierr = ISEqual(isrow,iscol,&isequal);CHKERRQ(ierr); 227 if (!isequal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices"); 228 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 229 230 ierr = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(ierr); 231 ierr = ISCreateStride(PetscObjectComm((PetscObject)mat),high-low,low,1,&ISrow);CHKERRQ(ierr); 232 ierr = MatCreateSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);CHKERRQ(ierr); 233 ierr = ISDestroy(&ISrow);CHKERRQ(ierr); 234 235 if (ctx->D1){ 236 ierr=VecDuplicate(ctx->D1,&D1);CHKERRQ(ierr); 237 ierr=VecCopy(ctx->D1,D1);CHKERRQ(ierr); 238 } else { 239 D1 = NULL; 240 } 241 242 if (ctx->D2){ 243 Vec D2sub; 244 245 ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);CHKERRQ(ierr); 246 ierr=VecDuplicate(D2sub,&D2);CHKERRQ(ierr); 247 ierr=VecCopy(D2sub,D2);CHKERRQ(ierr); 248 ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);CHKERRQ(ierr); 249 } else { 250 D2 = NULL; 251 } 252 253 ierr = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(ierr); 254 ierr = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(ierr); 255 ierr = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(ierr); 256 if (ctx->D1){ 257 ierr = PetscObjectDereference((PetscObject)D1);CHKERRQ(ierr); 258 } 259 if (ctx->D2){ 260 ierr = PetscObjectDereference((PetscObject)D2);CHKERRQ(ierr); 261 } 262 PetscFunctionReturn(0); 263 } 264 265 static PetscErrorCode MatCreateSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B) 266 { 267 PetscErrorCode ierr; 268 PetscInt i; 269 270 PetscFunctionBegin; 271 if (scall == MAT_INITIAL_MATRIX) { 272 ierr = PetscCalloc1(n+1,B );CHKERRQ(ierr); 273 } 274 for ( i=0; i<n; i++ ) { 275 ierr = MatCreateSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 276 } 277 PetscFunctionReturn(0); 278 } 279 280 static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col) 281 { 282 PetscErrorCode ierr; 283 PetscInt low,high; 284 PetscScalar zero=0.0,one=1.0; 285 286 PetscFunctionBegin; 287 ierr=VecSet(Y, zero);CHKERRQ(ierr); 288 ierr=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(ierr); 289 if (col>=low && col<high){ 290 ierr=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(ierr); 291 } 292 ierr=VecAssemblyBegin(Y);CHKERRQ(ierr); 293 ierr=VecAssemblyEnd(Y);CHKERRQ(ierr); 294 ierr=MatMult_ADA(mat,Y,Y);CHKERRQ(ierr); 295 PetscFunctionReturn(0); 296 } 297 298 PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat) 299 { 300 PetscErrorCode ierr; 301 PetscMPIInt size; 302 PetscBool sametype, issame, isdense, isseqdense; 303 TaoMatADACtx ctx; 304 305 PetscFunctionBegin; 306 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 307 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 308 309 ierr = PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr); 310 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);CHKERRQ(ierr); 311 ierr = PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);CHKERRQ(ierr); 312 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 313 314 if (sametype || issame) { 315 ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(ierr); 316 } else if (isdense) { 317 PetscInt i,j,low,high,m,n,M,N; 318 const PetscScalar *dptr; 319 Vec X; 320 321 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 322 ierr=MatGetSize(mat,&M,&N);CHKERRQ(ierr); 323 ierr=MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 324 ierr = MatCreateDense(PetscObjectComm((PetscObject)mat),m,m,N,N,NULL,NewMat);CHKERRQ(ierr); 325 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 326 for (i=0;i<M;i++){ 327 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 328 ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr); 329 for (j=0; j<high-low; j++){ 330 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 331 } 332 ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr); 333 } 334 ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 335 ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 336 ierr = VecDestroy(&X);CHKERRQ(ierr); 337 } else if (isseqdense && size==1){ 338 PetscInt i,j,low,high,m,n,M,N; 339 const PetscScalar *dptr; 340 Vec X; 341 342 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 343 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 344 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 345 ierr = MatCreateSeqDense(PetscObjectComm((PetscObject)mat),N,N,NULL,NewMat);CHKERRQ(ierr); 346 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 347 for (i=0;i<M;i++){ 348 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 349 ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr); 350 for (j=0; j<high-low; j++){ 351 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 352 } 353 ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr); 354 } 355 ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 356 ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 357 ierr = VecDestroy(&X);CHKERRQ(ierr); 358 } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type"); 359 PetscFunctionReturn(0); 360 } 361 362 static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm) 363 { 364 PetscErrorCode ierr; 365 TaoMatADACtx ctx; 366 367 PetscFunctionBegin; 368 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 369 if (type == NORM_FROBENIUS) { 370 *norm = 1.0; 371 } else if (type == NORM_1 || type == NORM_INFINITY) { 372 *norm = 1.0; 373 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 374 PetscFunctionReturn(0); 375 } 376 377 /*@C 378 MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal 379 380 Collective on matrix 381 382 Input Parameters: 383 + mat - matrix of arbitrary type 384 . d1 - A vector with diagonal elements of D1 385 - d2 - A vector 386 387 Output Parameters: 388 . J - New matrix whose operations are defined in terms of mat, D1, and D2. 389 390 Notes: 391 The user provides the input data and is responsible for destroying 392 this data after matrix J has been destroyed. 393 The operation MatMult(A,D2,D1) must be well defined. 394 Before calling the operation MatGetDiagonal(), the function 395 MatADAComputeDiagonal() must be called. The matrices A and D1 must 396 be the same during calls to MatADAComputeDiagonal() and 397 MatGetDiagonal(). 398 399 Level: developer 400 401 .seealso: MatCreate() 402 @*/ 403 PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J) 404 { 405 MPI_Comm comm=PetscObjectComm((PetscObject)mat); 406 TaoMatADACtx ctx; 407 PetscErrorCode ierr; 408 PetscInt nloc,n; 409 410 PetscFunctionBegin; 411 ierr = PetscNew(&ctx);CHKERRQ(ierr); 412 ctx->A=mat; 413 ctx->D1=d1; 414 ctx->D2=d2; 415 if (d1){ 416 ierr = VecDuplicate(d1,&ctx->W);CHKERRQ(ierr); 417 ierr = PetscObjectReference((PetscObject)d1);CHKERRQ(ierr); 418 } else { 419 ctx->W = NULL; 420 } 421 if (d2){ 422 ierr = VecDuplicate(d2,&ctx->W2);CHKERRQ(ierr); 423 ierr = VecDuplicate(d2,&ctx->ADADiag);CHKERRQ(ierr); 424 ierr = PetscObjectReference((PetscObject)d2);CHKERRQ(ierr); 425 } else { 426 ctx->W2 = NULL; 427 ctx->ADADiag = NULL; 428 } 429 430 ctx->GotDiag=0; 431 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 432 433 ierr=VecGetLocalSize(d2,&nloc);CHKERRQ(ierr); 434 ierr=VecGetSize(d2,&n);CHKERRQ(ierr); 435 436 ierr = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(ierr); 437 438 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);CHKERRQ(ierr); 439 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);CHKERRQ(ierr); 440 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);CHKERRQ(ierr); 441 ierr = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);CHKERRQ(ierr); 442 ierr = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);CHKERRQ(ierr); 443 ierr = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);CHKERRQ(ierr); 444 ierr = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);CHKERRQ(ierr); 445 ierr = MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);CHKERRQ(ierr); 446 ierr = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);CHKERRQ(ierr); 447 ierr = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);CHKERRQ(ierr); 448 ierr = MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_ADA);CHKERRQ(ierr); 449 ierr = MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);CHKERRQ(ierr); 450 ierr = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);CHKERRQ(ierr); 451 ierr = MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_ADA);CHKERRQ(ierr); 452 453 ierr = PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);CHKERRQ(ierr); 454 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));CHKERRQ(ierr); 455 456 ierr = MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 457 PetscFunctionReturn(0); 458 } 459