xref: /petsc/src/mat/tests/ex102.c (revision b122ec5aa1bd4469eb4e0673542fb7de3f411254)
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   PetscBool      flg;
18 
19   CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help));
20   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&M,NULL));
21   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
22 
23   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
24          Create the sparse matrix
25      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
26   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A));
27   CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N));
28   CHKERRQ(MatSetOptionsPrefix(A,"A_"));
29   CHKERRQ(MatSetFromOptions(A));
30   CHKERRQ(MatSeqAIJSetPreallocation(A,5,NULL));
31   CHKERRQ(MatMPIAIJSetPreallocation(A,5,NULL,5,NULL));
32   CHKERRQ(MatSetUp(A));
33   CHKERRQ(MatSetRandom(A,NULL));
34 
35   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36          Create the dense matrices
37      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&U));
39   CHKERRQ(MatSetSizes(U,PETSC_DECIDE,PETSC_DECIDE,M,3));
40   CHKERRQ(MatSetType(U,MATDENSE));
41   CHKERRQ(MatSetOptionsPrefix(U,"U_"));
42   CHKERRQ(MatSetFromOptions(U));
43   CHKERRQ(MatSetUp(U));
44   CHKERRQ(MatSetRandom(U,NULL));
45 
46   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&V));
47   CHKERRQ(MatSetSizes(V,PETSC_DECIDE,PETSC_DECIDE,N,3));
48   CHKERRQ(MatSetType(V,MATDENSE));
49   CHKERRQ(MatSetOptionsPrefix(V,"V_"));
50   CHKERRQ(MatSetFromOptions(V));
51   CHKERRQ(MatSetUp(V));
52   CHKERRQ(MatSetRandom(V,NULL));
53 
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   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-use_c",&flg));
61   if (flg) {
62     CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,3,&c));
63     CHKERRQ(VecSetRandom(c,NULL));
64   }
65 
66   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67          Create low rank correction matrix
68      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-low_rank",&flg));
70   if (flg) {
71     /* create a low-rank matrix, with no A-matrix */
72     CHKERRQ(MatCreateLRC(NULL,U,c,V,&LR));
73     CHKERRQ(MatDestroy(&A));
74   } else {
75     CHKERRQ(MatCreateLRC(A,U,c,V,&LR));
76   }
77 
78   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79          Create the low rank correction matrix explicitly to check for
80          correctness
81      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82   CHKERRQ(MatHermitianTranspose(V,MAT_INITIAL_MATRIX,&X));
83   CHKERRQ(MatDiagonalScale(X,c,NULL));
84   CHKERRQ(MatMatMult(U,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&LRe));
85   CHKERRQ(MatDestroy(&X));
86   if (A) {
87     CHKERRQ(MatAYPX(LRe,1.0,A,DIFFERENT_NONZERO_PATTERN));
88   }
89 
90   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91          Create test vectors
92      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93   CHKERRQ(MatCreateVecs(LR,&x,&b));
94   CHKERRQ(VecSetRandom(x,NULL));
95   CHKERRQ(MatMult(LR,x,b));
96   CHKERRQ(MatMultTranspose(LR,b,x));
97   CHKERRQ(VecDestroy(&x));
98   CHKERRQ(VecDestroy(&b));
99 
100   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101          Check correctness
102      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103   CHKERRQ(MatMultEqual(LR,LRe,10,&flg));
104   if (!flg) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMult\n"));
105 #if !defined(PETSC_USE_COMPLEX)
106   CHKERRQ(MatMultHermitianTransposeEqual(LR,LRe,10,&flg));
107   if (!flg) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMultTranspose\n"));
108 #endif
109 
110   CHKERRQ(MatDestroy(&A));
111   CHKERRQ(MatDestroy(&LRe));
112   CHKERRQ(MatDestroy(&U));
113   CHKERRQ(MatDestroy(&V));
114   CHKERRQ(VecDestroy(&c));
115   CHKERRQ(MatDestroy(&LR));
116 
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   CHKERRQ(PetscFinalize());
124   return 0;
125 }
126 
127 /*TEST
128 
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
139 
140 TEST*/
141