xref: /petsc/src/mat/tests/ex102.c (revision 40badf4fbc550ac1f60bd080eaff6de6d55b946d)
1 
2 static char help[] = "Tests MatCreateLRC()\n\n";
3 
4 /*T
5    Concepts: Low rank correction
6 
7    Processors: n
8 T*/
9 
10 #include <petscmat.h>
11 
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   PetscErrorCode ierr;
18   PetscBool      flg;
19 
20   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
21   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&M,NULL));
22   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
23 
24   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
25          Create the sparse matrix
26      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
27   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A));
28   CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N));
29   CHKERRQ(MatSetOptionsPrefix(A,"A_"));
30   CHKERRQ(MatSetFromOptions(A));
31   CHKERRQ(MatSeqAIJSetPreallocation(A,5,NULL));
32   CHKERRQ(MatMPIAIJSetPreallocation(A,5,NULL,5,NULL));
33   CHKERRQ(MatSetUp(A));
34   CHKERRQ(MatSetRandom(A,NULL));
35 
36   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37          Create the dense matrices
38      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&U));
40   CHKERRQ(MatSetSizes(U,PETSC_DECIDE,PETSC_DECIDE,M,3));
41   CHKERRQ(MatSetType(U,MATDENSE));
42   CHKERRQ(MatSetOptionsPrefix(U,"U_"));
43   CHKERRQ(MatSetFromOptions(U));
44   CHKERRQ(MatSetUp(U));
45   CHKERRQ(MatSetRandom(U,NULL));
46 
47   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&V));
48   CHKERRQ(MatSetSizes(V,PETSC_DECIDE,PETSC_DECIDE,N,3));
49   CHKERRQ(MatSetType(V,MATDENSE));
50   CHKERRQ(MatSetOptionsPrefix(V,"V_"));
51   CHKERRQ(MatSetFromOptions(V));
52   CHKERRQ(MatSetUp(V));
53   CHKERRQ(MatSetRandom(V,NULL));
54 
55   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56          Create a vector to hold the diagonal of C
57          A sequential vector can be created as well on each process
58          It is user responsibility to ensure the data in the vector
59          is consistent across processors
60      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-use_c",&flg));
62   if (flg) {
63     CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,3,&c));
64     CHKERRQ(VecSetRandom(c,NULL));
65   }
66 
67   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68          Create low rank correction matrix
69      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-low_rank",&flg));
71   if (flg) {
72     /* create a low-rank matrix, with no A-matrix */
73     CHKERRQ(MatCreateLRC(NULL,U,c,V,&LR));
74     CHKERRQ(MatDestroy(&A));
75   } else {
76     CHKERRQ(MatCreateLRC(A,U,c,V,&LR));
77   }
78 
79   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80          Create the low rank correction matrix explicitly to check for
81          correctness
82      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83   CHKERRQ(MatHermitianTranspose(V,MAT_INITIAL_MATRIX,&X));
84   CHKERRQ(MatDiagonalScale(X,c,NULL));
85   CHKERRQ(MatMatMult(U,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&LRe));
86   CHKERRQ(MatDestroy(&X));
87   if (A) {
88     CHKERRQ(MatAYPX(LRe,1.0,A,DIFFERENT_NONZERO_PATTERN));
89   }
90 
91   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92          Create test vectors
93      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94   CHKERRQ(MatCreateVecs(LR,&x,&b));
95   CHKERRQ(VecSetRandom(x,NULL));
96   CHKERRQ(MatMult(LR,x,b));
97   CHKERRQ(MatMultTranspose(LR,b,x));
98   CHKERRQ(VecDestroy(&x));
99   CHKERRQ(VecDestroy(&b));
100 
101   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102          Check correctness
103      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104   CHKERRQ(MatMultEqual(LR,LRe,10,&flg));
105   if (!flg) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMult\n"));
106 #if !defined(PETSC_USE_COMPLEX)
107   CHKERRQ(MatMultHermitianTransposeEqual(LR,LRe,10,&flg));
108   if (!flg) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMultTranspose\n"));
109 #endif
110 
111   CHKERRQ(MatDestroy(&A));
112   CHKERRQ(MatDestroy(&LRe));
113   CHKERRQ(MatDestroy(&U));
114   CHKERRQ(MatDestroy(&V));
115   CHKERRQ(VecDestroy(&c));
116   CHKERRQ(MatDestroy(&LR));
117 
118   /*
119      Always call PetscFinalize() before exiting a program.  This routine
120        - finalizes the PETSc libraries as well as MPI
121        - provides summary and diagnostic information if certain runtime
122          options are chosen (e.g., -log_view).
123   */
124   ierr = PetscFinalize();
125   return ierr;
126 }
127 
128 /*TEST
129 
130    testset:
131       output_file: output/ex102_1.out
132       nsize: {{1 2}}
133       args: -low_rank {{0 1}} -use_c {{0 1}}
134       test:
135         suffix: standard
136       test:
137         suffix: cuda
138         requires: cuda
139         args: -A_mat_type aijcusparse -U_mat_type densecuda -V_mat_type densecuda
140 
141 TEST*/
142