Actual source code: ex102.c
petsc-3.9.4 2018-09-11
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;
16: PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,rstart,rend;
18: PetscBool flg;
19: PetscScalar *u,a;
21: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
22: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
23: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
25: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
26: Create the sparse matrix
27: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
28: MatCreate(PETSC_COMM_WORLD,&A);
29: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
30: MatSetFromOptions(A);
31: MatSetUp(A);
32: MatGetOwnershipRange(A,&Istart,&Iend);
33: for (Ii=Istart; Ii<Iend; Ii++) {
34: a = -1.0; i = Ii/n; j = Ii - i*n;
35: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
36: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
37: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
38: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
39: a = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&a,INSERT_VALUES);
40: }
41: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
42: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
44: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: Create the dense matrices
46: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47: MatCreate(PETSC_COMM_WORLD,&U);
48: MatSetSizes(U,PETSC_DECIDE,PETSC_DECIDE,m*n,3);
49: MatSetType(U,MATDENSE);
50: MatSetUp(U);
51: MatGetOwnershipRange(U,&rstart,&rend);
52: MatDenseGetArray(U,&u);
53: for (i=rstart; i<rend; i++) {
54: u[i-rstart] = (PetscReal)i;
55: u[i+rend-2*rstart] = (PetscReal)1000*i;
56: u[i+2*rend-3*rstart] = (PetscReal)100000*i;
57: }
58: MatDenseRestoreArray(U,&u);
59: MatAssemblyBegin(U,MAT_FINAL_ASSEMBLY);
60: MatAssemblyEnd(U,MAT_FINAL_ASSEMBLY);
62: MatCreate(PETSC_COMM_WORLD,&V);
63: MatSetSizes(V,PETSC_DECIDE,PETSC_DECIDE,m*n,3);
64: MatSetType(V,MATDENSE);
65: MatSetUp(V);
66: MatGetOwnershipRange(U,&rstart,&rend);
67: MatDenseGetArray(V,&u);
68: for (i=rstart; i<rend; i++) {
69: u[i-rstart] = (PetscReal)i;
70: u[i+rend-2*rstart] = (PetscReal)1.2*i;
71: u[i+2*rend-3*rstart] = (PetscReal)1.67*i+2;
72: }
73: MatDenseRestoreArray(V,&u);
74: MatAssemblyBegin(V,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(V,MAT_FINAL_ASSEMBLY);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Create a vector to hold the diagonal of C
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: PetscOptionsHasName(NULL,NULL,"-use_c",&flg);
81: if (flg) {
82: VecCreateSeq(PETSC_COMM_SELF,3,&c);
83: VecGetArray(c,&u);
84: u[0] = 2.0;
85: u[1] = -1.0;
86: u[2] = 1.0;
87: VecRestoreArray(c,&u);
88: }
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Create low rank correction matrix
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: PetscOptionsHasName(NULL,NULL,"-low_rank",&flg);
94: if (flg) {
95: /* create a low-rank matrix, with no A-matrix */
96: MatCreateLRC(NULL,U,c,V,&LR);
97: } else {
98: MatCreateLRC(A,U,c,V,&LR);
99: }
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Create test vectors
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: VecCreate(PETSC_COMM_WORLD,&x);
105: VecSetSizes(x,PETSC_DECIDE,m*n);
106: VecSetFromOptions(x);
107: VecDuplicate(x,&b);
108: VecGetOwnershipRange(x,&rstart,&rend);
109: VecGetArray(x,&u);
110: for (i=rstart; i<rend; i++) u[i-rstart] = (PetscScalar)i;
111: VecRestoreArray(x,&u);
113: MatMult(LR,x,b);
114: /*
115: View the product if desired
116: */
117: PetscOptionsHasName(NULL,NULL,"-view_product",&flg);
118: if (flg) {VecView(b,PETSC_VIEWER_STDOUT_WORLD);}
120: VecDestroy(&x);
121: VecDestroy(&b);
122: /* you can destroy the matrices in any order you like */
123: MatDestroy(&A);
124: MatDestroy(&U);
125: MatDestroy(&V);
126: VecDestroy(&c);
127: MatDestroy(&LR);
129: /*
130: Always call PetscFinalize() before exiting a program. This routine
131: - finalizes the PETSc libraries as well as MPI
132: - provides summary and diagnostic information if certain runtime
133: options are chosen (e.g., -log_view).
134: */
135: PetscFinalize();
136: return ierr;
137: }
140: /*TEST
142: test:
143: suffix: 1
144: nsize: 2
145: args: -view_product
147: test:
148: suffix: 2
149: nsize: 2
150: args: -low_rank -view_product
152: test:
153: suffix: 3
154: nsize: 2
155: args: -use_c -view_product
157: TEST*/