Actual source code: ex32.c

slepc-3.15.1 2021-05-28
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 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";

 16: #include <slepclme.h>

 18: int main(int argc,char **argv)
 19: {
 20:   Mat                A;           /* problem matrix */
 21:   Mat                C,C1;        /* right-hand side */
 22:   Mat                X,X1;        /* solution */
 23:   LME                lme;
 24:   PetscReal          tol,errest,error;
 25:   PetscScalar        *u,sigma=0.0;
 26:   PetscInt           N,n=10,m,Istart,Iend,II,maxit,its,ncv,i,j,rank=0;
 27:   PetscErrorCode     ierr;
 28:   PetscBool          flag;
 29:   LMEConvergedReason reason;

 31:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 33:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 34:   PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag);
 35:   if (!flag) m=n;
 36:   N = n*m;
 37:   PetscOptionsGetScalar(NULL,NULL,"-sigma",&sigma,NULL);
 38:   PetscOptionsGetInt(NULL,NULL,"-rank",&rank,NULL);
 39:   PetscPrintf(PETSC_COMM_WORLD,"\nLyapunov equation, N=%D (%Dx%D grid)\n\n",N,n,m);

 41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 42:                        Create the 2-D Laplacian, A
 43:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 45:   MatCreate(PETSC_COMM_WORLD,&A);
 46:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
 47:   MatSetFromOptions(A);
 48:   MatSetUp(A);
 49:   MatGetOwnershipRange(A,&Istart,&Iend);
 50:   for (II=Istart;II<Iend;II++) {
 51:     i = II/n; j = II-i*n;
 52:     if (i>0) { MatSetValue(A,II,II-n,1.0,INSERT_VALUES); }
 53:     if (i<m-1) { MatSetValue(A,II,II+n,1.0,INSERT_VALUES); }
 54:     if (j>0) { MatSetValue(A,II,II-1,1.0,INSERT_VALUES); }
 55:     if (j<n-1) { MatSetValue(A,II,II+1,1.0,INSERT_VALUES); }
 56:     MatSetValue(A,II,II,-4.0-sigma,INSERT_VALUES);
 57:   }
 58:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 59:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 61:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62:        Create a low-rank Mat to store the right-hand side C = C1*C1'
 63:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 65:   MatCreate(PETSC_COMM_WORLD,&C1);
 66:   MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,N,2);
 67:   MatSetType(C1,MATDENSE);
 68:   MatSetUp(C1);
 69:   MatGetOwnershipRange(C1,&Istart,&Iend);
 70:   MatDenseGetArray(C1,&u);
 71:   for (i=Istart;i<Iend;i++) {
 72:     if (i<N/2) u[i-Istart] = 1.0;
 73:     if (i==0) u[i+Iend-2*Istart] = -2.0;
 74:     if (i==1) u[i+Iend-2*Istart] = -1.0;
 75:     if (i==2) u[i+Iend-2*Istart] = -1.0;
 76:   }
 77:   MatDenseRestoreArray(C1,&u);
 78:   MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY);
 79:   MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);
 80:   MatCreateLRC(NULL,C1,NULL,NULL,&C);
 81:   MatDestroy(&C1);

 83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84:                 Create the solver and set various options
 85:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 86:   /*
 87:      Create the matrix equation solver context
 88:   */
 89:   LMECreate(PETSC_COMM_WORLD,&lme);

 91:   /*
 92:      Set the type of equation
 93:   */
 94:   LMESetProblemType(lme,LME_LYAPUNOV);

 96:   /*
 97:      Set the matrix coefficients, the right-hand side, and the solution.
 98:      In this case, it is a Lyapunov equation A*X+X*A'=-C where both
 99:      C and X are symmetric and low-rank, C=C1*C1', X=X1*X1'
100:   */
101:   LMESetCoefficients(lme,A,NULL,NULL,NULL);
102:   LMESetRHS(lme,C);

104:   if (rank) {  /* Create X only if the user has specified a nonzero value of rank */
105:     PetscPrintf(PETSC_COMM_WORLD," Computing a solution with prescribed rank=%d\n",rank);
106:     MatCreate(PETSC_COMM_WORLD,&X1);
107:     MatSetSizes(X1,PETSC_DECIDE,PETSC_DECIDE,N,rank);
108:     MatSetType(X1,MATDENSE);
109:     MatSetUp(X1);
110:     MatAssemblyBegin(X1,MAT_FINAL_ASSEMBLY);
111:     MatAssemblyEnd(X1,MAT_FINAL_ASSEMBLY);
112:     MatCreateLRC(NULL,X1,NULL,NULL,&X);
113:     MatDestroy(&X1);
114:     LMESetSolution(lme,X);
115:     MatDestroy(&X);
116:   }

118:   /*
119:      (Optional) Set other solver options
120:   */
121:   LMESetTolerances(lme,1e-07,PETSC_DEFAULT);

123:   /*
124:      Set solver parameters at runtime
125:   */
126:   LMESetFromOptions(lme);

128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:                    Solve the matrix equation, A*X+X*A'=-C
130:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

132:   LMESolve(lme);
133:   LMEGetConvergedReason(lme,&reason);
134:   if (reason<0) SETERRQ(PETSC_COMM_WORLD,1,"Solver did not converge");

136:   if (!rank) {  /* X1 was created by the solver, so extract it and see how many columns it has */
137:     LMEGetSolution(lme,&X);
138:     MatLRCGetMats(X,NULL,&X1,NULL,NULL);
139:     MatGetSize(X1,NULL,&rank);
140:     PetscPrintf(PETSC_COMM_WORLD," The solver has computed a solution with rank=%d\n",rank);
141:   }

143:   /*
144:      Optional: Get some information from the solver and display it
145:   */
146:   LMEGetIterationNumber(lme,&its);
147:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
148:   LMEGetDimensions(lme,&ncv);
149:   PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %D\n",ncv);
150:   LMEGetTolerances(lme,&tol,&maxit);
151:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);

153:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154:                         Compute residual error
155:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

157:   LMEGetErrorEstimate(lme,&errest);
158:   PetscPrintf(PETSC_COMM_WORLD," Error estimate reported by the solver: %.4g\n",(double)errest);
159:   if (n<=150) {
160:     LMEComputeError(lme,&error);
161:     PetscPrintf(PETSC_COMM_WORLD," Computed residual norm: %.4g\n\n",(double)error);
162:   } else {
163:     PetscPrintf(PETSC_COMM_WORLD," Matrix too large to compute residual norm\n\n");
164:   }

166:   /*
167:      Free work space
168:   */
169:   LMEDestroy(&lme);
170:   MatDestroy(&A);
171:   MatDestroy(&C);
172:   SlepcFinalize();
173:   return ierr;
174: }

176: /*TEST

178:    test:
179:       suffix: 1
180:       requires: !single

182:    test:
183:       suffix: 2
184:       args: -rank 40
185:       requires: !single

187: TEST*/