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,Mat *B) 153 { 154 PetscErrorCode ierr; 155 TaoMatADACtx ctx; 156 157 PetscFunctionBegin; 158 if (*B){ 159 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 160 ierr = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(ierr); 161 } 162 PetscFunctionReturn(0); 163 } 164 165 PetscErrorCode MatADAComputeDiagonal(Mat mat) 166 { 167 PetscErrorCode ierr; 168 PetscInt i,m,n,low,high; 169 PetscScalar *dtemp,*dptr; 170 TaoMatADACtx ctx; 171 172 PetscFunctionBegin; 173 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 174 ierr = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(ierr); 175 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 176 177 ierr = PetscMalloc1(n,&dtemp );CHKERRQ(ierr); 178 179 for (i=0; i<n; i++){ 180 ierr = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(ierr); 181 ierr = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(ierr); 182 ierr = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 183 } 184 for (i=0; i<n; i++){ 185 ierr = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr); 186 } 187 188 ierr = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 189 for (i=low; i<high; i++){ 190 dptr[i-low]= dtemp[i]; 191 } 192 ierr = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(ierr); 193 ierr = PetscFree(dtemp);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 197 static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v) 198 { 199 PetscErrorCode ierr; 200 PetscReal one=1.0; 201 TaoMatADACtx ctx; 202 203 PetscFunctionBegin; 204 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 205 ierr = MatADAComputeDiagonal(mat);CHKERRQ(ierr); 206 ierr=VecCopy(ctx->ADADiag,v);CHKERRQ(ierr); 207 if (ctx->D2){ 208 ierr=VecAXPY(v, one, ctx->D2);CHKERRQ(ierr); 209 } 210 PetscFunctionReturn(0); 211 } 212 213 static PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat) 214 { 215 PetscErrorCode ierr; 216 PetscInt low,high; 217 IS ISrow; 218 Vec D1,D2; 219 Mat Atemp; 220 TaoMatADACtx ctx; 221 PetscBool isequal; 222 223 PetscFunctionBegin; 224 ierr = ISEqual(isrow,iscol,&isequal);CHKERRQ(ierr); 225 if (!isequal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices"); 226 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 227 228 ierr = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(ierr); 229 ierr = ISCreateStride(PetscObjectComm((PetscObject)mat),high-low,low,1,&ISrow);CHKERRQ(ierr); 230 ierr = MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);CHKERRQ(ierr); 231 ierr = ISDestroy(&ISrow);CHKERRQ(ierr); 232 233 if (ctx->D1){ 234 ierr=VecDuplicate(ctx->D1,&D1);CHKERRQ(ierr); 235 ierr=VecCopy(ctx->D1,D1);CHKERRQ(ierr); 236 } else { 237 D1 = NULL; 238 } 239 240 if (ctx->D2){ 241 Vec D2sub; 242 243 ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);CHKERRQ(ierr); 244 ierr=VecDuplicate(D2sub,&D2);CHKERRQ(ierr); 245 ierr=VecCopy(D2sub,D2);CHKERRQ(ierr); 246 ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);CHKERRQ(ierr); 247 } else { 248 D2 = NULL; 249 } 250 251 ierr = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(ierr); 252 ierr = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(ierr); 253 ierr = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(ierr); 254 if (ctx->D1){ 255 ierr = PetscObjectDereference((PetscObject)D1);CHKERRQ(ierr); 256 } 257 if (ctx->D2){ 258 ierr = PetscObjectDereference((PetscObject)D2);CHKERRQ(ierr); 259 } 260 PetscFunctionReturn(0); 261 } 262 263 static PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B) 264 { 265 PetscErrorCode ierr; 266 PetscInt i; 267 268 PetscFunctionBegin; 269 if (scall == MAT_INITIAL_MATRIX) { 270 ierr = PetscMalloc1(n+1,B );CHKERRQ(ierr); 271 } 272 for ( i=0; i<n; i++ ) { 273 ierr = MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 274 } 275 PetscFunctionReturn(0); 276 } 277 278 static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col) 279 { 280 PetscErrorCode ierr; 281 PetscInt low,high; 282 PetscScalar zero=0.0,one=1.0; 283 284 PetscFunctionBegin; 285 ierr=VecSet(Y, zero);CHKERRQ(ierr); 286 ierr=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(ierr); 287 if (col>=low && col<high){ 288 ierr=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(ierr); 289 } 290 ierr=VecAssemblyBegin(Y);CHKERRQ(ierr); 291 ierr=VecAssemblyEnd(Y);CHKERRQ(ierr); 292 ierr=MatMult_ADA(mat,Y,Y);CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat) 297 { 298 PetscErrorCode ierr; 299 PetscMPIInt size; 300 PetscBool sametype, issame, isdense, isseqdense; 301 TaoMatADACtx ctx; 302 303 PetscFunctionBegin; 304 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 305 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 306 307 ierr = PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr); 308 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);CHKERRQ(ierr); 309 ierr = PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);CHKERRQ(ierr); 310 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 311 312 if (sametype || issame) { 313 ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(ierr); 314 } else if (isdense) { 315 PetscInt i,j,low,high,m,n,M,N; 316 const PetscScalar *dptr; 317 Vec X; 318 319 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 320 ierr=MatGetSize(mat,&M,&N);CHKERRQ(ierr); 321 ierr=MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 322 ierr = MatCreateDense(PetscObjectComm((PetscObject)mat),m,m,N,N,NULL,NewMat);CHKERRQ(ierr); 323 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 324 for (i=0;i<M;i++){ 325 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 326 ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr); 327 for (j=0; j<high-low; j++){ 328 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 329 } 330 ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr); 331 } 332 ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 333 ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 334 ierr = VecDestroy(&X);CHKERRQ(ierr); 335 } else if (isseqdense && size==1){ 336 PetscInt i,j,low,high,m,n,M,N; 337 const PetscScalar *dptr; 338 Vec X; 339 340 ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr); 341 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 342 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 343 ierr = MatCreateSeqDense(PetscObjectComm((PetscObject)mat),N,N,NULL,NewMat);CHKERRQ(ierr); 344 ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr); 345 for (i=0;i<M;i++){ 346 ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr); 347 ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr); 348 for (j=0; j<high-low; j++){ 349 ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr); 350 } 351 ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr); 352 } 353 ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 354 ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 355 ierr = VecDestroy(&X);CHKERRQ(ierr); 356 } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type"); 357 PetscFunctionReturn(0); 358 } 359 360 static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm) 361 { 362 PetscErrorCode ierr; 363 TaoMatADACtx ctx; 364 365 PetscFunctionBegin; 366 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 367 if (type == NORM_FROBENIUS) { 368 *norm = 1.0; 369 } else if (type == NORM_1 || type == NORM_INFINITY) { 370 *norm = 1.0; 371 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 372 PetscFunctionReturn(0); 373 } 374 375 /*@C 376 MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal 377 378 Collective on matrix 379 380 Input Parameters: 381 + mat - matrix of arbitrary type 382 . d1 - A vector with diagonal elements of D1 383 - d2 - A vector 384 385 Output Parameters: 386 . J - New matrix whose operations are defined in terms of mat, D1, and D2. 387 388 Notes: 389 The user provides the input data and is responsible for destroying 390 this data after matrix J has been destroyed. 391 The operation MatMult(A,D2,D1) must be well defined. 392 Before calling the operation MatGetDiagonal(), the function 393 MatADAComputeDiagonal() must be called. The matrices A and D1 must 394 be the same during calls to MatADAComputeDiagonal() and 395 MatGetDiagonal(). 396 397 Level: developer 398 399 .seealso: MatCreate() 400 @*/ 401 PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J) 402 { 403 MPI_Comm comm=PetscObjectComm((PetscObject)mat); 404 TaoMatADACtx ctx; 405 PetscErrorCode ierr; 406 PetscInt nloc,n; 407 408 PetscFunctionBegin; 409 ierr = PetscNew(&ctx);CHKERRQ(ierr); 410 ctx->A=mat; 411 ctx->D1=d1; 412 ctx->D2=d2; 413 if (d1){ 414 ierr = VecDuplicate(d1,&ctx->W);CHKERRQ(ierr); 415 ierr = PetscObjectReference((PetscObject)d1);CHKERRQ(ierr); 416 } else { 417 ctx->W = NULL; 418 } 419 if (d2){ 420 ierr = VecDuplicate(d2,&ctx->W2);CHKERRQ(ierr); 421 ierr = VecDuplicate(d2,&ctx->ADADiag);CHKERRQ(ierr); 422 ierr = PetscObjectReference((PetscObject)d2);CHKERRQ(ierr); 423 } else { 424 ctx->W2 = NULL; 425 ctx->ADADiag = NULL; 426 } 427 428 ctx->GotDiag=0; 429 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 430 431 ierr=VecGetLocalSize(d2,&nloc);CHKERRQ(ierr); 432 ierr=VecGetSize(d2,&n);CHKERRQ(ierr); 433 434 ierr = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(ierr); 435 436 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);CHKERRQ(ierr); 437 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);CHKERRQ(ierr); 438 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);CHKERRQ(ierr); 439 ierr = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);CHKERRQ(ierr); 440 ierr = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);CHKERRQ(ierr); 441 ierr = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);CHKERRQ(ierr); 442 ierr = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);CHKERRQ(ierr); 443 ierr = MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);CHKERRQ(ierr); 444 ierr = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);CHKERRQ(ierr); 445 ierr = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);CHKERRQ(ierr); 446 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);CHKERRQ(ierr); 447 ierr = MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);CHKERRQ(ierr); 448 ierr = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);CHKERRQ(ierr); 449 ierr = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);CHKERRQ(ierr); 450 451 ierr = PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);CHKERRQ(ierr); 452 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));CHKERRQ(ierr); 453 454 ierr = MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 455 PetscFunctionReturn(0); 456 } 457