Actual source code: ex102.c
petsc-3.3-p6 2013-02-11
2: /* Program usage: mpiexec -n <procs> ex2 [-help] [all PETSc options] */
4: static char help[] = "Tests MatCreateLRC()\n\n";
6: /*T
7: Concepts: Low rank correction
9: Processors: n
10: T*/
12: #include <petscmat.h>
16: int main(int argc,char **args)
17: {
18: Vec x,b;
19: Mat A,U,V,LR;
20: PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,rstart,rend;
22: PetscBool flg;
23: PetscScalar *u,a;
25: PetscInitialize(&argc,&args,(char *)0,help);
26: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
27: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
29: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: Create the sparse matrix
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
32: MatCreate(PETSC_COMM_WORLD,&A);
33: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
34: MatSetFromOptions(A);
35: MatSetUp(A);
36: MatGetOwnershipRange(A,&Istart,&Iend);
37: for (Ii=Istart; Ii<Iend; Ii++) {
38: a = -1.0; i = Ii/n; j = Ii - i*n;
39: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
40: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
41: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
42: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);}
43: a = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&a,INSERT_VALUES);
44: }
45: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
46: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Create the dense matrices
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: MatCreate(PETSC_COMM_WORLD,&U);
52: MatSetSizes(U,PETSC_DECIDE,PETSC_DECIDE,m*n,3);
53: MatSetType(U,MATDENSE);
54: MatSetUp(U);
55: MatGetOwnershipRange(U,&rstart,&rend);
56: MatGetArray(U,&u);
57: for (i=rstart; i<rend; i++) {
58: u[i-rstart] = (PetscReal)i;
59: u[i+rend-2*rstart] = (PetscReal)1000*i;
60: u[i+2*rend-3*rstart] = (PetscReal)100000*i;
61: }
62: MatRestoreArray(U,&u);
63: MatAssemblyBegin(U,MAT_FINAL_ASSEMBLY);
64: MatAssemblyEnd(U,MAT_FINAL_ASSEMBLY);
67: MatCreate(PETSC_COMM_WORLD,&V);
68: MatSetSizes(V,PETSC_DECIDE,PETSC_DECIDE,m*n,3);
69: MatSetType(V,MATDENSE);
70: MatSetUp(V);
71: MatGetOwnershipRange(U,&rstart,&rend);
72: MatGetArray(V,&u);
73: for (i=rstart; i<rend; i++) {
74: u[i-rstart] = (PetscReal)i;
75: u[i+rend-2*rstart] = (PetscReal)1.2*i;
76: u[i+2*rend-3*rstart] = (PetscReal)1.67*i+2;
77: }
78: MatRestoreArray(V,&u);
79: MatAssemblyBegin(V,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(V,MAT_FINAL_ASSEMBLY);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Create low rank created matrix
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: MatCreateLRC(A,U,V,&LR);
85: MatSetUp(LR);
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Create test vectors
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: VecCreate(PETSC_COMM_WORLD,&x);
91: VecSetSizes(x,PETSC_DECIDE,m*n);
92: VecSetFromOptions(x);
93: VecDuplicate(x,&b);
94: VecGetOwnershipRange(x,&rstart,&rend);
95: VecGetArray(x,&u);
96: for (i=rstart; i<rend; i++) u[i-rstart] = (PetscScalar)i;
97: VecRestoreArray(x,&u);
99: MatMult(LR,x,b);
100: /*
101: View the product if desired
102: */
103: PetscOptionsHasName(PETSC_NULL,"-view_product",&flg);
104: if (flg) {VecView(b,PETSC_VIEWER_STDOUT_WORLD);}
106: VecDestroy(&x);
107: VecDestroy(&b);
108: /* you can destroy the matrices in any order you like */
109: MatDestroy(&A);
110: MatDestroy(&U);
111: MatDestroy(&V);
112: MatDestroy(&LR);
114: /*
115: Always call PetscFinalize() before exiting a program. This routine
116: - finalizes the PETSc libraries as well as MPI
117: - provides summary and diagnostic information if certain runtime
118: options are chosen (e.g., -log_summary).
119: */
120: PetscFinalize();
121: return 0;
122: }