Actual source code: ex26.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: static char help[] = "Tests MatGetRowIJ for SeqAIJ, SeqBAIJ and SeqSBAIJ\n\n";

  4: #include <petscmat.h>

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

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

 39: int main(int argc,char **args)
 40: {
 41:   Mat            A,B,C;
 42:   PetscInt       i,j,k,m = 3,n = 3,bs = 1;
 44:   PetscMPIInt    size;

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

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

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

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