1 #include <petscmat.h> /*I "mat.h" I*/ 2 3 PETSC_INTERN PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*); 4 5 typedef struct { 6 Mat A; 7 Vec D1; 8 Vec D2; 9 Vec W; 10 Vec W2; 11 Vec ADADiag; 12 PetscInt GotDiag; 13 } _p_TaoMatADACtx; 14 typedef _p_TaoMatADACtx* TaoMatADACtx; 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,&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(Mat M,Vec D, InsertMode mode) 46 { 47 TaoMatADACtx ctx; 48 PetscReal zero=0.0,one = 1.0; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)M),PETSC_ERR_SUP,"Cannot insert diagonal entries of this matrix type, can only add"); 53 ierr = MatShellGetContext(M,&ctx);CHKERRQ(ierr); 54 if (!ctx->D2) { 55 ierr = VecDuplicate(D,&ctx->D2);CHKERRQ(ierr); 56 ierr = VecSet(ctx->D2, zero);CHKERRQ(ierr); 57 } 58 ierr = VecAXPY(ctx->D2, one, D);CHKERRQ(ierr); 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode MatDestroy_ADA(Mat mat) 63 { 64 PetscErrorCode ierr; 65 TaoMatADACtx ctx; 66 67 PetscFunctionBegin; 68 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 69 ierr = VecDestroy(&ctx->W);CHKERRQ(ierr); 70 ierr = VecDestroy(&ctx->W2);CHKERRQ(ierr); 71 ierr = VecDestroy(&ctx->ADADiag);CHKERRQ(ierr); 72 ierr = MatDestroy(&ctx->A);CHKERRQ(ierr); 73 ierr = VecDestroy(&ctx->D1);CHKERRQ(ierr); 74 ierr = VecDestroy(&ctx->D2);CHKERRQ(ierr); 75 ierr = PetscFree(ctx);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer) 80 { 81 PetscFunctionBegin; 82 PetscFunctionReturn(0); 83 } 84 85 static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a) 86 { 87 PetscErrorCode ierr; 88 TaoMatADACtx ctx; 89 90 PetscFunctionBegin; 91 ierr = MatShellGetContext(Y,&ctx);CHKERRQ(ierr); 92 ierr = VecShift(ctx->D2,a);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M) 97 { 98 PetscErrorCode ierr; 99 TaoMatADACtx ctx; 100 Mat A2; 101 Vec D1b=NULL,D2b; 102 103 PetscFunctionBegin; 104 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 105 ierr = MatDuplicate(ctx->A,op,&A2);CHKERRQ(ierr); 106 if (ctx->D1) { 107 ierr = VecDuplicate(ctx->D1,&D1b);CHKERRQ(ierr); 108 ierr = VecCopy(ctx->D1,D1b);CHKERRQ(ierr); 109 } 110 ierr = VecDuplicate(ctx->D2,&D2b);CHKERRQ(ierr); 111 ierr = VecCopy(ctx->D2,D2b);CHKERRQ(ierr); 112 ierr = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(ierr); 113 if (ctx->D1) { 114 ierr = PetscObjectDereference((PetscObject)D1b);CHKERRQ(ierr); 115 } 116 ierr = PetscObjectDereference((PetscObject)D2b);CHKERRQ(ierr); 117 ierr = PetscObjectDereference((PetscObject)A2);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg) 122 { 123 PetscErrorCode ierr; 124 TaoMatADACtx ctx1,ctx2; 125 126 PetscFunctionBegin; 127 ierr = MatShellGetContext(A,&ctx1);CHKERRQ(ierr); 128 ierr = MatShellGetContext(B,&ctx2);CHKERRQ(ierr); 129 ierr = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(ierr); 130 if (*flg==PETSC_TRUE) { 131 ierr = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(ierr); 132 } 133 if (*flg==PETSC_TRUE) { 134 ierr = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(ierr); 135 } 136 PetscFunctionReturn(0); 137 } 138 139 static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a) 140 { 141 PetscErrorCode ierr; 142 TaoMatADACtx ctx; 143 144 PetscFunctionBegin; 145 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 146 ierr = VecScale(ctx->D1,a);CHKERRQ(ierr); 147 if (ctx->D2) { 148 ierr = VecScale(ctx->D2,a);CHKERRQ(ierr); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 static PetscErrorCode MatTranspose_ADA(Mat mat,MatReuse reuse,Mat *B) 154 { 155 PetscErrorCode ierr; 156 TaoMatADACtx ctx; 157 158 PetscFunctionBegin; 159 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 160 if (reuse == MAT_INITIAL_MATRIX) { 161 ierr = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(ierr); 162 } else if (reuse == MAT_REUSE_MATRIX) { 163 ierr = MatCopy(mat,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 164 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Does not support inplace transpose"); 165 PetscFunctionReturn(0); 166 } 167 168 static PetscErrorCode MatADAComputeDiagonal(Mat mat) 169 { 170 PetscErrorCode ierr; 171 PetscInt i,m,n,low,high; 172 PetscScalar *dtemp,*dptr; 173 TaoMatADACtx ctx; 174 175 PetscFunctionBegin; 176 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 177 ierr = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(ierr); 178 ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr); 179 180 ierr = PetscMalloc1(n,&dtemp);CHKERRQ(ierr); 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,&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 identical column and row indices"); 228 ierr = MatShellGetContext(mat,&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,&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,&ctx);CHKERRQ(ierr); 307 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRMPI(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,PETSC_ERR_PLIB,"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,&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 defining a diagonal matrix 385 - d2 - A vector defining a diagonal matrix 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 394 Level: developer 395 396 .seealso: MatCreate() 397 @*/ 398 PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J) 399 { 400 MPI_Comm comm = PetscObjectComm((PetscObject)mat); 401 TaoMatADACtx ctx; 402 PetscErrorCode ierr; 403 PetscInt nloc,n; 404 405 PetscFunctionBegin; 406 ierr = PetscNew(&ctx);CHKERRQ(ierr); 407 ctx->A=mat; 408 ctx->D1=d1; 409 ctx->D2=d2; 410 if (d1) { 411 ierr = VecDuplicate(d1,&ctx->W);CHKERRQ(ierr); 412 ierr = PetscObjectReference((PetscObject)d1);CHKERRQ(ierr); 413 } else { 414 ctx->W = NULL; 415 } 416 if (d2) { 417 ierr = VecDuplicate(d2,&ctx->W2);CHKERRQ(ierr); 418 ierr = VecDuplicate(d2,&ctx->ADADiag);CHKERRQ(ierr); 419 ierr = PetscObjectReference((PetscObject)d2);CHKERRQ(ierr); 420 } else { 421 ctx->W2 = NULL; 422 ctx->ADADiag = NULL; 423 } 424 425 ctx->GotDiag = 0; 426 ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr); 427 428 ierr = VecGetLocalSize(d2,&nloc);CHKERRQ(ierr); 429 ierr = VecGetSize(d2,&n);CHKERRQ(ierr); 430 431 ierr = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(ierr); 432 ierr = MatShellSetManageScalingShifts(*J);CHKERRQ(ierr); 433 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);CHKERRQ(ierr); 434 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);CHKERRQ(ierr); 435 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);CHKERRQ(ierr); 436 ierr = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);CHKERRQ(ierr); 437 ierr = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);CHKERRQ(ierr); 438 ierr = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);CHKERRQ(ierr); 439 ierr = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);CHKERRQ(ierr); 440 ierr = MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);CHKERRQ(ierr); 441 ierr = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);CHKERRQ(ierr); 442 ierr = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);CHKERRQ(ierr); 443 ierr = MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_ADA);CHKERRQ(ierr); 444 ierr = MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);CHKERRQ(ierr); 445 ierr = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);CHKERRQ(ierr); 446 ierr = MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_ADA);CHKERRQ(ierr); 447 448 ierr = PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);CHKERRQ(ierr); 449 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));CHKERRQ(ierr); 450 451 ierr = MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454