Actual source code: ex102.c


  2: static char help[] = "Tests MatCreateLRC()\n\n";

  4: /*T
  5:    Concepts: Low rank correction

  7:    Processors: n
  8: T*/

 10: #include <petscmat.h>

 12: int main(int argc,char **args)
 13: {
 14:   Vec            x,b,c=NULL;
 15:   Mat            A,U,V,LR,X,LRe;
 16:   PetscInt       M = 5, N = 7;
 17:   PetscBool      flg;

 19:   PetscInitialize(&argc,&args,(char*)0,help);
 20:   PetscOptionsGetInt(NULL,NULL,"-m",&M,NULL);
 21:   PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);

 23:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 24:          Create the sparse matrix
 25:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 26:   MatCreate(PETSC_COMM_WORLD,&A);
 27:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
 28:   MatSetOptionsPrefix(A,"A_");
 29:   MatSetFromOptions(A);
 30:   MatSeqAIJSetPreallocation(A,5,NULL);
 31:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
 32:   MatSetUp(A);
 33:   MatSetRandom(A,NULL);

 35:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 36:          Create the dense matrices
 37:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 38:   MatCreate(PETSC_COMM_WORLD,&U);
 39:   MatSetSizes(U,PETSC_DECIDE,PETSC_DECIDE,M,3);
 40:   MatSetType(U,MATDENSE);
 41:   MatSetOptionsPrefix(U,"U_");
 42:   MatSetFromOptions(U);
 43:   MatSetUp(U);
 44:   MatSetRandom(U,NULL);

 46:   MatCreate(PETSC_COMM_WORLD,&V);
 47:   MatSetSizes(V,PETSC_DECIDE,PETSC_DECIDE,N,3);
 48:   MatSetType(V,MATDENSE);
 49:   MatSetOptionsPrefix(V,"V_");
 50:   MatSetFromOptions(V);
 51:   MatSetUp(V);
 52:   MatSetRandom(V,NULL);

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:          Create a vector to hold the diagonal of C
 56:          A sequential vector can be created as well on each process
 57:          It is user responsibility to ensure the data in the vector
 58:          is consistent across processors
 59:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 60:   PetscOptionsHasName(NULL,NULL,"-use_c",&flg);
 61:   if (flg) {
 62:     VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,3,&c);
 63:     VecSetRandom(c,NULL);
 64:   }

 66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:          Create low rank correction matrix
 68:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 69:   PetscOptionsHasName(NULL,NULL,"-low_rank",&flg);
 70:   if (flg) {
 71:     /* create a low-rank matrix, with no A-matrix */
 72:     MatCreateLRC(NULL,U,c,V,&LR);
 73:     MatDestroy(&A);
 74:   } else {
 75:     MatCreateLRC(A,U,c,V,&LR);
 76:   }

 78:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79:          Create the low rank correction matrix explicitly to check for
 80:          correctness
 81:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 82:   MatHermitianTranspose(V,MAT_INITIAL_MATRIX,&X);
 83:   MatDiagonalScale(X,c,NULL);
 84:   MatMatMult(U,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&LRe);
 85:   MatDestroy(&X);
 86:   if (A) {
 87:     MatAYPX(LRe,1.0,A,DIFFERENT_NONZERO_PATTERN);
 88:   }

 90:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 91:          Create test vectors
 92:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 93:   MatCreateVecs(LR,&x,&b);
 94:   VecSetRandom(x,NULL);
 95:   MatMult(LR,x,b);
 96:   MatMultTranspose(LR,b,x);
 97:   VecDestroy(&x);
 98:   VecDestroy(&b);

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:          Check correctness
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   MatMultEqual(LR,LRe,10,&flg);
104:   if (!flg) PetscPrintf(PETSC_COMM_WORLD,"Error in MatMult\n");
105: #if !defined(PETSC_USE_COMPLEX)
106:   MatMultHermitianTransposeEqual(LR,LRe,10,&flg);
107:   if (!flg) PetscPrintf(PETSC_COMM_WORLD,"Error in MatMultTranspose\n");
108: #endif

110:   MatDestroy(&A);
111:   MatDestroy(&LRe);
112:   MatDestroy(&U);
113:   MatDestroy(&V);
114:   VecDestroy(&c);
115:   MatDestroy(&LR);

117:   /*
118:      Always call PetscFinalize() before exiting a program.  This routine
119:        - finalizes the PETSc libraries as well as MPI
120:        - provides summary and diagnostic information if certain runtime
121:          options are chosen (e.g., -log_view).
122:   */
123:   PetscFinalize();
124:   return 0;
125: }

127: /*TEST

129:    testset:
130:       output_file: output/ex102_1.out
131:       nsize: {{1 2}}
132:       args: -low_rank {{0 1}} -use_c {{0 1}}
133:       test:
134:         suffix: standard
135:       test:
136:         suffix: cuda
137:         requires: cuda
138:         args: -A_mat_type aijcusparse -U_mat_type densecuda -V_mat_type densecuda

140: TEST*/