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