Actual source code: ex2.c

      1: /*
      2:  * Example code testing SeqDense matrices with an LDA (leading dimension
      3:  * of the user-allocated arrray) larger than M.
      4:  * This example tests the functionality of MatInsertValues, MatMult,
      5:  * and MatMultTranspose.
      6:  */
      7: #include <stdlib.h>
     8:  #include petscmat.h

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

     21:   PetscInitialize(&argc,&argv,0,0);

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

     40:   PetscMalloc(size*sizeof(PetscScalar),&x);
     41:   for (i=0; i<size; i++) {
     42:     x[i] = one;
     43:   }
     44:   VecCreateSeqWithArray(MPI_COMM_SELF,size,x,&X);
     45:   VecAssemblyBegin(X);
     46:   VecAssemblyEnd(X);

     48:   PetscMalloc(size*sizeof(PetscScalar),&y);
     49:   VecCreateSeqWithArray(MPI_COMM_SELF,size,y,&Y);
     50:   VecAssemblyBegin(Y);
     51:   VecAssemblyEnd(Y);

     53:   PetscMalloc(size*sizeof(PetscScalar),&z);
     54:   VecCreateSeqWithArray(MPI_COMM_SELF,size,z,&Z);
     55:   VecAssemblyBegin(Z);
     56:   VecAssemblyEnd(Z);

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

     69:   MatCreate(MPI_COMM_SELF,&A12);
     70:   MatSetSizes(A12,size1,size2,size1,size2);
     71:   MatSetType(A12,MATSEQDENSE);
     72:   MatSeqDenseSetPreallocation(A12,b+size1*size);
     73:   MatSeqDenseSetLDA(A12,size);
     74:   MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
     75:   MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);

     77:   MatCreate(MPI_COMM_SELF,&A21);
     78:   MatSetSizes(A21,size2,size1,size2,size1);
     79:   MatSetType(A21,MATSEQDENSE);
     80:   MatSeqDenseSetPreallocation(A21,b+size1);
     81:   MatSeqDenseSetLDA(A21,size);
     82:   MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
     83:   MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);

     85:   MatCreate(MPI_COMM_SELF,&A22);
     86:   MatSetSizes(A22,size2,size2,size2,size2);
     87:   MatSetType(A22,MATSEQDENSE);
     88:   MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
     89:   MatSeqDenseSetLDA(A22,size);
     90:   MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
     91:   MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);

     93:   VecCreateSeqWithArray(MPI_COMM_SELF,size1,x,&X1);
     94:   VecCreateSeqWithArray(MPI_COMM_SELF,size2,x+size1,&X2);
     95:   VecCreateSeqWithArray(MPI_COMM_SELF,size1,z,&Z1);
     96:   VecCreateSeqWithArray(MPI_COMM_SELF,size2,z+size1,&Z2);

     98:   /*
     99:    * Now multiple matrix times input in two ways;
    100:    * compare the results
    101:    */
    102:   MatMult(A,X,Y);
    103:   MatMult(A11,X1,Z1);
    104:   MatMultAdd(A12,X2,Z1,Z1);
    105:   MatMult(A22,X2,Z2);
    106:   MatMultAdd(A21,X1,Z2,Z2);
    107:   VecAXPY(Z,mone,Y);
    108:   VecNorm(Z,NORM_2,&nrm);
    109:   printf("Test1; error norm=%e\n",nrm);
    110:   /*
    111:   printf("MatMult the usual way:\n"); VecView(Y,0);
    112:   printf("MatMult by subblock:\n"); VecView(Z,0);
    113:   */

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

    133:   MatMult(A,X,Y);
    134:   MatMult(A11,X1,Z1);
    135:   MatMultAdd(A12,X2,Z1,Z1);
    136:   MatMult(A22,X2,Z2);
    137:   MatMultAdd(A21,X1,Z2,Z2);
    138:   VecAXPY(Z,mone,Y);
    139:   VecNorm(Z,NORM_2,&nrm);
    140:   printf("Test2; error norm=%e\n",nrm);

    142:   /*
    143:    * Transpose product
    144:    */
    145:   MatMultTranspose(A,X,Y);
    146:   MatMultTranspose(A11,X1,Z1);
    147:   MatMultTransposeAdd(A21,X2,Z1,Z1);
    148:   MatMultTranspose(A22,X2,Z2);
    149:   MatMultTransposeAdd(A12,X1,Z2,Z2);
    150:   VecAXPY(Z,mone,Y);
    151:   VecNorm(Z,NORM_2,&nrm);
    152:   printf("Test3; error norm=%e\n",nrm);

    154:   PetscFinalize();
    155:   return 0;
    156: }