Actual source code: ex102.c


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

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Vec       x, b, c = NULL;
  9:   Mat       A, U, V, LR, X, LRe;
 10:   PetscInt  M = 5, N = 7;
 11:   PetscBool flg;

 14:   PetscInitialize(&argc, &args, (char *)0, help);
 15:   PetscOptionsGetInt(NULL, NULL, "-m", &M, NULL);
 16:   PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL);

 18:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 19:          Create the sparse matrix
 20:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 21:   MatCreate(PETSC_COMM_WORLD, &A);
 22:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
 23:   MatSetOptionsPrefix(A, "A_");
 24:   MatSetFromOptions(A);
 25:   MatSeqAIJSetPreallocation(A, 5, NULL);
 26:   MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
 27:   MatSetUp(A);
 28:   MatSetRandom(A, NULL);

 30:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 31:          Create the dense matrices
 32:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 33:   MatCreate(PETSC_COMM_WORLD, &U);
 34:   MatSetSizes(U, PETSC_DECIDE, PETSC_DECIDE, M, 3);
 35:   MatSetType(U, MATDENSE);
 36:   MatSetOptionsPrefix(U, "U_");
 37:   MatSetFromOptions(U);
 38:   MatSetUp(U);
 39:   MatSetRandom(U, NULL);

 41:   MatCreate(PETSC_COMM_WORLD, &V);
 42:   MatSetSizes(V, PETSC_DECIDE, PETSC_DECIDE, N, 3);
 43:   MatSetType(V, MATDENSE);
 44:   MatSetOptionsPrefix(V, "V_");
 45:   MatSetFromOptions(V);
 46:   MatSetUp(V);
 47:   MatSetRandom(V, NULL);

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

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

 73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:          Create the low rank correction matrix explicitly to check for
 75:          correctness
 76:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 77:   MatHermitianTranspose(V, MAT_INITIAL_MATRIX, &X);
 78:   MatDiagonalScale(X, c, NULL);
 79:   MatMatMult(U, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &LRe);
 80:   MatDestroy(&X);
 81:   if (A) MatAYPX(LRe, 1.0, A, DIFFERENT_NONZERO_PATTERN);

 83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84:          Create test vectors
 85:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 86:   MatCreateVecs(LR, &x, &b);
 87:   VecSetRandom(x, NULL);
 88:   MatMult(LR, x, b);
 89:   MatMultTranspose(LR, b, x);
 90:   VecDestroy(&x);
 91:   VecDestroy(&b);

 93:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:          Check correctness
 95:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 96:   MatMultEqual(LR, LRe, 10, &flg);
 97:   if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error in MatMult\n");
 98: #if !defined(PETSC_USE_COMPLEX)
 99:   MatMultHermitianTransposeEqual(LR, LRe, 10, &flg);
100:   if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error in MatMultTranspose\n");
101: #endif

103:   MatDestroy(&A);
104:   MatDestroy(&LRe);
105:   MatDestroy(&U);
106:   MatDestroy(&V);
107:   VecDestroy(&c);
108:   MatDestroy(&LR);

110:   /*
111:      Always call PetscFinalize() before exiting a program.  This routine
112:        - finalizes the PETSc libraries as well as MPI
113:        - provides summary and diagnostic information if certain runtime
114:          options are chosen (e.g., -log_view).
115:   */
116:   PetscFinalize();
117:   return 0;
118: }

120: /*TEST

122:    testset:
123:       output_file: output/ex102_1.out
124:       nsize: {{1 2}}
125:       args: -low_rank {{0 1}} -use_c {{0 1}}
126:       test:
127:         suffix: standard
128:       test:
129:         suffix: cuda
130:         requires: cuda
131:         args: -A_mat_type aijcusparse -U_mat_type densecuda -V_mat_type densecuda

133: TEST*/