Actual source code: ex97.c
petsc-3.4.2 2013-07-02

      1: static const char help[] = "Tests MatGetSubMatrix with MatSubMatrix versus MatAIJ, non-square\n";

      3: #include <petscmat.h>

      7: static PetscErrorCode AssembleMatrix(MPI_Comm comm,Mat *A)
      8: {
     10:   Mat            B;
     11:   PetscInt       i,ms,me;

     14:   MatCreate(comm,&B);
     15:   MatSetSizes(B,5,6,PETSC_DETERMINE,PETSC_DETERMINE);
     16:   MatSetFromOptions(B);
     17:   MatSetUp(B);
     18:   MatGetOwnershipRange(B,&ms,&me);
     19:   for (i=ms; i<me; i++) {
     20:     MatSetValue(B,i,i,1.0*i,INSERT_VALUES);
     21:   }
     22:   MatSetValue(B,me-1,me,me*me,INSERT_VALUES);
     23:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
     24:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
     25:   *A   = B;
     26:   return(0);
     27: }

     31: static PetscErrorCode Compare2(Vec *X,const char *test)
     32: {
     34:   PetscReal      norm;
     35:   Vec            Y;
     36:   PetscInt       verbose = 0;

     39:   VecDuplicate(X[0],&Y);
     40:   VecCopy(X[0],Y);
     41:   VecAYPX(Y,-1.0,X[1]);
     42:   VecNorm(Y,NORM_INFINITY,&norm);

     44:   PetscOptionsGetInt(NULL,"-verbose",&verbose,NULL);
     45:   if (norm < 1.e-12 && verbose < 1) {
     46:     PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference < 1e-12\n",test);
     47:   } else {
     48:     PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference %G\n",test,norm);
     49:   }
     50:   if (verbose > 1) {
     51:     VecView(X[0],PETSC_VIEWER_STDOUT_WORLD);
     52:     VecView(X[1],PETSC_VIEWER_STDOUT_WORLD);
     53:     VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
     54:   }
     55:   VecDestroy(&Y);
     56:   return(0);
     57: }

     61: static PetscErrorCode CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1)
     62: {
     64:   Vec            *ltmp,*rtmp;

     67:   VecDuplicateVecs(right,2,&rtmp);
     68:   VecDuplicateVecs(left,2,&ltmp);
     69:   MatScale(A,PETSC_PI);
     70:   MatScale(B,PETSC_PI);
     71:   MatDiagonalScale(A,left,right);
     72:   MatDiagonalScale(B,left,right);

     74:   MatMult(A,X,ltmp[0]);
     75:   MatMult(B,X,ltmp[1]);
     76:   Compare2(ltmp,"MatMult");

     78:   MatMultTranspose(A,Y,rtmp[0]);
     79:   MatMultTranspose(B,Y,rtmp[1]);
     80:   Compare2(rtmp,"MatMultTranspose");

     82:   VecCopy(Y1,ltmp[0]);
     83:   VecCopy(Y1,ltmp[1]);
     84:   MatMultAdd(A,X,ltmp[0],ltmp[0]);
     85:   MatMultAdd(B,X,ltmp[1],ltmp[1]);
     86:   Compare2(ltmp,"MatMultAdd v2==v3");

     88:   MatMultAdd(A,X,Y1,ltmp[0]);
     89:   MatMultAdd(B,X,Y1,ltmp[1]);
     90:   Compare2(ltmp,"MatMultAdd v2!=v3");

     92:   VecCopy(X1,rtmp[0]);
     93:   VecCopy(X1,rtmp[1]);
     94:   MatMultTransposeAdd(A,Y,rtmp[0],rtmp[0]);
     95:   MatMultTransposeAdd(B,Y,rtmp[1],rtmp[1]);
     96:   Compare2(rtmp,"MatMultTransposeAdd v2==v3");

     98:   MatMultTransposeAdd(A,Y,X1,rtmp[0]);
     99:   MatMultTransposeAdd(B,Y,X1,rtmp[1]);
    100:   Compare2(rtmp,"MatMultTransposeAdd v2!=v3");

    102:   VecDestroyVecs(2,&ltmp);
    103:   VecDestroyVecs(2,&rtmp);
    104:   return(0);
    105: }

    109: int main(int argc, char *argv[])
    110: {
    112:   Mat            A,B,Asub,Bsub;
    113:   PetscInt       ms,idxrow[3],idxcol[4];
    114:   Vec            left,right,X,Y,X1,Y1;
    115:   IS             isrow,iscol;
    116:   PetscBool      random = PETSC_TRUE;

    118:   PetscInitialize(&argc,&argv,NULL,help);
    119:   AssembleMatrix(PETSC_COMM_WORLD,&A);
    120:   AssembleMatrix(PETSC_COMM_WORLD,&B);
    121:   MatShellSetOperation(B,MATOP_GET_SUBMATRIX,NULL);
    122:   MatShellSetOperation(B,MATOP_GET_SUBMATRICES,NULL);
    123:   MatGetOwnershipRange(A,&ms,NULL);

    125:   idxrow[0] = ms+1;
    126:   idxrow[1] = ms+2;
    127:   idxrow[2] = ms+4;
    128:   ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);

    130:   idxcol[0] = ms+1;
    131:   idxcol[1] = ms+2;
    132:   idxcol[2] = ms+4;
    133:   idxcol[3] = ms+5;
    134:   ISCreateGeneral(PETSC_COMM_WORLD,4,idxcol,PETSC_USE_POINTER,&iscol);

    136:   MatGetSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);
    137:   MatGetSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);

    139:   MatGetVecs(Asub,&right,&left);
    140:   VecDuplicate(right,&X);
    141:   VecDuplicate(right,&X1);
    142:   VecDuplicate(left,&Y);
    143:   VecDuplicate(left,&Y1);

    145:   PetscOptionsGetBool(NULL,"-random",&random,NULL);
    146:   if (random) {
    147:     VecSetRandom(right,NULL);
    148:     VecSetRandom(left,NULL);
    149:     VecSetRandom(X,NULL);
    150:     VecSetRandom(Y,NULL);
    151:     VecSetRandom(X1,NULL);
    152:     VecSetRandom(Y1,NULL);
    153:   } else {
    154:     VecSet(right,1.0);
    155:     VecSet(left,2.0);
    156:     VecSet(X,3.0);
    157:     VecSet(Y,4.0);
    158:     VecSet(X1,3.0);
    159:     VecSet(Y1,4.0);
    160:   }
    161:   CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1);
    162:   ISDestroy(&isrow);
    163:   ISDestroy(&iscol);
    164:   MatDestroy(&A);
    165:   MatDestroy(&B);
    166:   MatDestroy(&Asub);
    167:   MatDestroy(&Bsub);
    168:   VecDestroy(&left);
    169:   VecDestroy(&right);
    170:   VecDestroy(&X);
    171:   VecDestroy(&Y);
    172:   VecDestroy(&X1);
    173:   VecDestroy(&Y1);
    174:   PetscFinalize();
    175:   return 0;
    176: }