Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Solves a Lypunov equation with the shifted 2-D Laplacian.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
14 : " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
15 :
16 : #include <slepclme.h>
17 :
18 2 : int main(int argc,char **argv)
19 : {
20 2 : Mat A; /* problem matrix */
21 2 : Mat C,C1; /* right-hand side */
22 2 : Mat X,X1; /* solution */
23 2 : LME lme;
24 2 : PetscReal tol,errest,error;
25 2 : PetscScalar *u,sigma=0.0;
26 2 : PetscInt N,n=10,m,Istart,Iend,II,maxit,its,ncv,i,j,rank=0;
27 2 : PetscBool flag;
28 2 : LMEConvergedReason reason;
29 :
30 2 : PetscFunctionBeginUser;
31 2 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
32 :
33 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
34 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
35 2 : if (!flag) m=n;
36 2 : N = n*m;
37 2 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-sigma",&sigma,NULL));
38 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-rank",&rank,NULL));
39 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nLyapunov equation, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
40 :
41 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42 : Create the 2-D Laplacian, A
43 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44 :
45 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
46 2 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
47 2 : PetscCall(MatSetFromOptions(A));
48 2 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
49 202 : for (II=Istart;II<Iend;II++) {
50 200 : i = II/n; j = II-i*n;
51 200 : if (i>0) PetscCall(MatSetValue(A,II,II-n,1.0,INSERT_VALUES));
52 200 : if (i<m-1) PetscCall(MatSetValue(A,II,II+n,1.0,INSERT_VALUES));
53 200 : if (j>0) PetscCall(MatSetValue(A,II,II-1,1.0,INSERT_VALUES));
54 200 : if (j<n-1) PetscCall(MatSetValue(A,II,II+1,1.0,INSERT_VALUES));
55 200 : PetscCall(MatSetValue(A,II,II,-4.0-sigma,INSERT_VALUES));
56 : }
57 2 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
58 2 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
59 :
60 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61 : Create a low-rank Mat to store the right-hand side C = C1*C1'
62 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63 :
64 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C1));
65 2 : PetscCall(MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,N,2));
66 2 : PetscCall(MatSetType(C1,MATDENSE));
67 2 : PetscCall(MatGetOwnershipRange(C1,&Istart,&Iend));
68 2 : PetscCall(MatDenseGetArray(C1,&u));
69 202 : for (i=Istart;i<Iend;i++) {
70 200 : if (i<N/2) u[i-Istart] = 1.0;
71 200 : if (i==0) u[i+Iend-2*Istart] = -2.0;
72 200 : if (i==1) u[i+Iend-2*Istart] = -1.0;
73 200 : if (i==2) u[i+Iend-2*Istart] = -1.0;
74 : }
75 2 : PetscCall(MatDenseRestoreArray(C1,&u));
76 2 : PetscCall(MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY));
77 2 : PetscCall(MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY));
78 2 : PetscCall(MatCreateLRC(NULL,C1,NULL,NULL,&C));
79 2 : PetscCall(MatDestroy(&C1));
80 :
81 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 : Create the solver and set various options
83 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84 : /*
85 : Create the matrix equation solver context
86 : */
87 2 : PetscCall(LMECreate(PETSC_COMM_WORLD,&lme));
88 :
89 : /*
90 : Set the type of equation
91 : */
92 2 : PetscCall(LMESetProblemType(lme,LME_LYAPUNOV));
93 :
94 : /*
95 : Set the matrix coefficients, the right-hand side, and the solution.
96 : In this case, it is a Lyapunov equation A*X+X*A'=-C where both
97 : C and X are symmetric and low-rank, C=C1*C1', X=X1*X1'
98 : */
99 2 : PetscCall(LMESetCoefficients(lme,A,NULL,NULL,NULL));
100 2 : PetscCall(LMESetRHS(lme,C));
101 :
102 2 : if (rank) { /* Create X only if the user has specified a nonzero value of rank */
103 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computing a solution with prescribed rank=%" PetscInt_FMT "\n",rank));
104 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&X1));
105 1 : PetscCall(MatSetSizes(X1,PETSC_DECIDE,PETSC_DECIDE,N,rank));
106 1 : PetscCall(MatSetType(X1,MATDENSE));
107 1 : PetscCall(MatAssemblyBegin(X1,MAT_FINAL_ASSEMBLY));
108 1 : PetscCall(MatAssemblyEnd(X1,MAT_FINAL_ASSEMBLY));
109 1 : PetscCall(MatCreateLRC(NULL,X1,NULL,NULL,&X));
110 1 : PetscCall(MatDestroy(&X1));
111 1 : PetscCall(LMESetSolution(lme,X));
112 1 : PetscCall(MatDestroy(&X));
113 : }
114 :
115 : /*
116 : (Optional) Set other solver options
117 : */
118 2 : PetscCall(LMESetTolerances(lme,1e-07,PETSC_CURRENT));
119 :
120 : /*
121 : Set solver parameters at runtime
122 : */
123 2 : PetscCall(LMESetFromOptions(lme));
124 :
125 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126 : Solve the matrix equation, A*X+X*A'=-C
127 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128 :
129 2 : PetscCall(LMESolve(lme));
130 2 : PetscCall(LMEGetConvergedReason(lme,&reason));
131 2 : PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
132 :
133 2 : if (!rank) { /* X1 was created by the solver, so extract it and see how many columns it has */
134 1 : PetscCall(LMEGetSolution(lme,&X));
135 1 : PetscCall(MatLRCGetMats(X,NULL,&X1,NULL,NULL));
136 1 : PetscCall(MatGetSize(X1,NULL,&rank));
137 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," The solver has computed a solution with rank=%" PetscInt_FMT "\n",rank));
138 : }
139 :
140 : /*
141 : Optional: Get some information from the solver and display it
142 : */
143 2 : PetscCall(LMEGetIterationNumber(lme,&its));
144 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its));
145 2 : PetscCall(LMEGetDimensions(lme,&ncv));
146 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv));
147 2 : PetscCall(LMEGetTolerances(lme,&tol,&maxit));
148 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
149 :
150 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151 : Compute residual error
152 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153 :
154 2 : PetscCall(LMEGetErrorEstimate(lme,&errest));
155 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Error estimate reported by the solver: %.4g\n",(double)errest));
156 2 : if (n<=150) {
157 2 : PetscCall(LMEComputeError(lme,&error));
158 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed residual norm: %.4g\n\n",(double)error));
159 0 : } else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix too large to compute residual norm\n\n"));
160 :
161 : /*
162 : Free work space
163 : */
164 2 : PetscCall(LMEDestroy(&lme));
165 2 : PetscCall(MatDestroy(&A));
166 2 : PetscCall(MatDestroy(&C));
167 2 : PetscCall(SlepcFinalize());
168 : return 0;
169 : }
170 :
171 : /*TEST
172 :
173 : test:
174 : suffix: 1
175 : requires: !single
176 :
177 : test:
178 : suffix: 2
179 : args: -rank 40
180 : requires: !single
181 :
182 : TEST*/
|
