Actual source code: ex26.c

petsc-3.8.4 2018-03-24
Report Typos and Errors

  2: static char help[] = "Tests MatGetRowIJ for SeqAIJ, SeqBAIJ and SeqSBAIJ\n\n";

  4:  #include <petscmat.h>

  6: PetscErrorCode DumpCSR(Mat A,PetscInt shift,PetscBool symmetric,PetscBool compressed)
  7: {
  8:   MatType        type;
  9:   PetscInt       i,j,nr,bs = 1;
 10:   const PetscInt *ia,*ja;
 11:   PetscBool      done,isseqbaij,isseqsbaij;

 15:   PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isseqbaij);
 16:   PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isseqsbaij);
 17:   if (isseqbaij || isseqsbaij) {
 18:     MatGetBlockSize(A,&bs);
 19:   }
 20:   MatGetType(A,&type);
 21:   MatGetRowIJ(A,shift,symmetric,compressed,&nr,&ia,&ja,&done);
 22:   PetscPrintf(PETSC_COMM_SELF,"===========================================================\n");
 23:   PetscPrintf(PETSC_COMM_SELF,"CSR for %s: shift %D symmetric %D compressed %D\n",type,shift,(PetscInt)symmetric,(PetscInt)compressed);
 24:   for (i=0;i<nr;i++) {
 25:     PetscPrintf(PETSC_COMM_SELF,"%D:",i+shift);
 26:     for (j=ia[i];j<ia[i+1];j++) {
 27:       PetscPrintf(PETSC_COMM_SELF," %D",ja[j-shift]);
 28:     }
 29:     PetscPrintf(PETSC_COMM_SELF,"\n");
 30:   }
 31:   MatRestoreRowIJ(A,shift,symmetric,compressed,&nr,&ia,&ja,&done);
 32:   return(0);
 33: }

 35: int main(int argc,char **args)
 36: {
 37:   Mat            A,B,C;
 38:   PetscInt       i,j,k,m = 3,n = 3,bs = 1;
 40:   PetscMPIInt    size;

 42:   PetscInitialize(&argc,&args,(char*)0,help);
 43:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 44:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 45:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 46:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 47:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 48:   /* adjust sizes by block size */
 49:   if (m%bs) m += bs-m%bs;
 50:   if (n%bs) n += bs-n%bs;

 52:   MatCreate(PETSC_COMM_SELF,&A);
 53:   MatSetSizes(A,m*n,m*n,PETSC_DECIDE,PETSC_DECIDE);
 54:   MatSetBlockSize(A,bs);
 55:   MatSetType(A,MATSEQAIJ);
 56:   MatSetUp(A);
 57:   MatCreate(PETSC_COMM_SELF,&B);
 58:   MatSetSizes(B,m*n,m*n,PETSC_DECIDE,PETSC_DECIDE);
 59:   MatSetBlockSize(B,bs);
 60:   MatSetType(B,MATSEQBAIJ);
 61:   MatSetUp(B);
 62:   MatCreate(PETSC_COMM_SELF,&C);
 63:   MatSetSizes(C,m*n,m*n,PETSC_DECIDE,PETSC_DECIDE);
 64:   MatSetBlockSize(C,bs);
 65:   MatSetType(C,MATSEQSBAIJ);
 66:   MatSetUp(C);
 67:   MatSetOption(C,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);

 69:   for (i=0; i<m; i++) {
 70:     for (j=0; j<n; j++) {
 71: 
 72:       PetscScalar v = -1.0;
 73:       PetscInt    Ii = j + n*i,J;
 74:       J = Ii - n;
 75:       if (J>=0)  {
 76:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 77:         MatSetValues(B,1,&Ii,1,&J,&v,INSERT_VALUES);
 78:         MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
 79:       }
 80:       J = Ii + n;
 81:       if (J<m*n) {
 82:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 83:         MatSetValues(B,1,&Ii,1,&J,&v,INSERT_VALUES);
 84:         MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
 85:       }
 86:       J = Ii - 1;
 87:       if (J>=0)  {
 88:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 89:         MatSetValues(B,1,&Ii,1,&J,&v,INSERT_VALUES);
 90:         MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
 91:       }
 92:       J = Ii + 1;
 93:       if (J<m*n) {
 94:         MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 95:         MatSetValues(B,1,&Ii,1,&J,&v,INSERT_VALUES);
 96:         MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
 97:       }
 98:       v = 4.0;
 99:       MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
100:       MatSetValues(B,1,&Ii,1,&Ii,&v,INSERT_VALUES);
101:       MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
102:     }
103:   }
104:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
105:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
106:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
107:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
108:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
109:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

111:   /* test MatGetRowIJ for the three Mat types */
112:   MatView(A,NULL);
113:   MatView(B,NULL);
114:   MatView(C,NULL);
115:   for (i=0;i<2;i++) {
116:     PetscInt shift = i;
117:     for (j=0;j<2;j++) {
118:       PetscBool symmetric = ((j>0) ? PETSC_FALSE : PETSC_TRUE);
119:       for (k=0;k<2;k++) {
120:         PetscBool compressed = ((k>0) ? PETSC_FALSE : PETSC_TRUE);
121:         DumpCSR(A,shift,symmetric,compressed);
122:         DumpCSR(B,shift,symmetric,compressed);
123:         DumpCSR(C,shift,symmetric,compressed);
124:       }
125:     }
126:   }
127:   MatDestroy(&A);
128:   MatDestroy(&B);
129:   MatDestroy(&C);
130:   PetscFinalize();
131:   return ierr;
132: }