Actual source code: ex2.c

      1: static char help[] = "testing SeqDense matrices with an LDA (leading dimension of the user-allocated arrray) larger than M.\n";
      2: /*
      3:  * Example code testing SeqDense matrices with an LDA (leading dimension
      4:  * of the user-allocated arrray) larger than M.
      5:  * This example tests the functionality of MatSetValues(), MatMult(),
      6:  * and MatMultTranspose().
      7:  */

      9: #include <petscmat.h>

     11: int main(int argc,char **argv)
     12: {
     13:   Mat            A,A11,A12,A21,A22;
     14:   Vec            X,X1,X2,Y,Z,Z1,Z2;
     15:   PetscScalar    *a,*b,*x,*y,*z,v,one=1;
     16:   PetscReal      nrm;
     18:   PetscInt       size=8,size1=6,size2=2, i,j;
     19:   PetscRandom    rnd;

     21:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
     22:   PetscRandomCreate(PETSC_COMM_SELF,&rnd);

     24:   /*
     25:    * Create matrix and three vectors: these are all normal
     26:    */
     27:   PetscMalloc1(size*size,&a);
     28:   PetscMalloc1(size*size,&b);
     29:   for (i=0; i<size; i++) {
     30:     for (j=0; j<size; j++) {
     31:       PetscRandomGetValue(rnd,&a[i+j*size]);
     32:       b[i+j*size] = a[i+j*size];
     33:     }
     34:   }
     35:   MatCreate(PETSC_COMM_SELF,&A);
     36:   MatSetSizes(A,size,size,size,size);
     37:   MatSetType(A,MATSEQDENSE);
     38:   MatSeqDenseSetPreallocation(A,a);
     39:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     40:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

     42:   PetscMalloc1(size,&x);
     43:   for (i=0; i<size; i++) {
     44:     x[i] = one;
     45:   }
     46:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,size,x,&X);
     47:   VecAssemblyBegin(X);
     48:   VecAssemblyEnd(X);

     50:   PetscMalloc1(size,&y);
     51:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,size,y,&Y);
     52:   VecAssemblyBegin(Y);
     53:   VecAssemblyEnd(Y);

     55:   PetscMalloc1(size,&z);
     56:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,size,z,&Z);
     57:   VecAssemblyBegin(Z);
     58:   VecAssemblyEnd(Z);

     60:   /*
     61:    * Now create submatrices and subvectors
     62:    */
     63:   MatCreate(PETSC_COMM_SELF,&A11);
     64:   MatSetSizes(A11,size1,size1,size1,size1);
     65:   MatSetType(A11,MATSEQDENSE);
     66:   MatSeqDenseSetPreallocation(A11,b);
     67:   MatDenseSetLDA(A11,size);
     68:   MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);
     69:   MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);

     71:   MatCreate(PETSC_COMM_SELF,&A12);
     72:   MatSetSizes(A12,size1,size2,size1,size2);
     73:   MatSetType(A12,MATSEQDENSE);
     74:   MatSeqDenseSetPreallocation(A12,b+size1*size);
     75:   MatDenseSetLDA(A12,size);
     76:   MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
     77:   MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);

     79:   MatCreate(PETSC_COMM_SELF,&A21);
     80:   MatSetSizes(A21,size2,size1,size2,size1);
     81:   MatSetType(A21,MATSEQDENSE);
     82:   MatSeqDenseSetPreallocation(A21,b+size1);
     83:   MatDenseSetLDA(A21,size);
     84:   MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
     85:   MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);

     87:   MatCreate(PETSC_COMM_SELF,&A22);
     88:   MatSetSizes(A22,size2,size2,size2,size2);
     89:   MatSetType(A22,MATSEQDENSE);
     90:   MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
     91:   MatDenseSetLDA(A22,size);
     92:   MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
     93:   MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);

     95:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,size1,x,&X1);
     96:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,size2,x+size1,&X2);
     97:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,size1,z,&Z1);
     98:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,size2,z+size1,&Z2);

    100:   /*
    101:    * Now multiple matrix times input in two ways;
    102:    * compare the results
    103:    */
    104:   MatMult(A,X,Y);
    105:   MatMult(A11,X1,Z1);
    106:   MatMultAdd(A12,X2,Z1,Z1);
    107:   MatMult(A22,X2,Z2);
    108:   MatMultAdd(A21,X1,Z2,Z2);
    109:   VecAXPY(Z,-1.0,Y);
    110:   VecNorm(Z,NORM_2,&nrm);
    111:   if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
    112:     PetscPrintf(PETSC_COMM_WORLD,"Test1; error norm=%g\n",(double)nrm);
    113:   }

    115:   /*
    116:    * Next test: change both matrices
    117:    */
    118:   PetscRandomGetValue(rnd,&v);
    119:   i    = 1;
    120:   j    = size-2;
    121:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
    122:   j   -= size1;
    123:   MatSetValues(A12,1,&i,1,&j,&v,INSERT_VALUES);
    124:   PetscRandomGetValue(rnd,&v);
    125:   i    = j = size1+1;
    126:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
    127:   i     =j = 1;
    128:   MatSetValues(A22,1,&i,1,&j,&v,INSERT_VALUES);
    129:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    130:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    131:   MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
    132:   MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
    133:   MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
    134:   MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);

    136:   MatMult(A,X,Y);
    137:   MatMult(A11,X1,Z1);
    138:   MatMultAdd(A12,X2,Z1,Z1);
    139:   MatMult(A22,X2,Z2);
    140:   MatMultAdd(A21,X1,Z2,Z2);
    141:   VecAXPY(Z,-1.0,Y);
    142:   VecNorm(Z,NORM_2,&nrm);
    143:   if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
    144:     PetscPrintf(PETSC_COMM_WORLD,"Test2; error norm=%g\n",(double)nrm);
    145:   }

    147:   /*
    148:    * Transpose product
    149:    */
    150:   MatMultTranspose(A,X,Y);
    151:   MatMultTranspose(A11,X1,Z1);
    152:   MatMultTransposeAdd(A21,X2,Z1,Z1);
    153:   MatMultTranspose(A22,X2,Z2);
    154:   MatMultTransposeAdd(A12,X1,Z2,Z2);
    155:   VecAXPY(Z,-1.0,Y);
    156:   VecNorm(Z,NORM_2,&nrm);
    157:   if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
    158:     PetscPrintf(PETSC_COMM_WORLD,"Test3; error norm=%g\n",(double)nrm);
    159:   }
    160:   PetscFree(a);
    161:   PetscFree(b);
    162:   PetscFree(x);
    163:   PetscFree(y);
    164:   PetscFree(z);
    165:   MatDestroy(&A);
    166:   MatDestroy(&A11);
    167:   MatDestroy(&A12);
    168:   MatDestroy(&A21);
    169:   MatDestroy(&A22);

    171:   VecDestroy(&X);
    172:   VecDestroy(&Y);
    173:   VecDestroy(&Z);

    175:   VecDestroy(&X1);
    176:   VecDestroy(&X2);
    177:   VecDestroy(&Z1);
    178:   VecDestroy(&Z2);
    179:   PetscRandomDestroy(&rnd);
    180:   PetscFinalize();
    181:   return ierr;
    182: }

    185: /*TEST

    187:    test:

    189: TEST*/