Actual source code: ex32.c

slepc-3.20.1 2023-11-27
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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:   PetscBool          flag;
 28:   LMEConvergedReason reason;

 30:   PetscFunctionBeginUser;
 31:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));

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

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

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

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

 65:   PetscCall(MatCreate(PETSC_COMM_WORLD,&C1));
 66:   PetscCall(MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,N,2));
 67:   PetscCall(MatSetType(C1,MATDENSE));
 68:   PetscCall(MatSetUp(C1));
 69:   PetscCall(MatGetOwnershipRange(C1,&Istart,&Iend));
 70:   PetscCall(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:   PetscCall(MatDenseRestoreArray(C1,&u));
 78:   PetscCall(MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY));
 79:   PetscCall(MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY));
 80:   PetscCall(MatCreateLRC(NULL,C1,NULL,NULL,&C));
 81:   PetscCall(MatDestroy(&C1));

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

 91:   /*
 92:      Set the type of equation
 93:   */
 94:   PetscCall(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:   PetscCall(LMESetCoefficients(lme,A,NULL,NULL,NULL));
102:   PetscCall(LMESetRHS(lme,C));

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

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

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

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

132:   PetscCall(LMESolve(lme));
133:   PetscCall(LMEGetConvergedReason(lme,&reason));
134:   PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");

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

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

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

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

164:   /*
165:      Free work space
166:   */
167:   PetscCall(LMEDestroy(&lme));
168:   PetscCall(MatDestroy(&A));
169:   PetscCall(MatDestroy(&C));
170:   PetscCall(SlepcFinalize());
171:   return 0;
172: }

174: /*TEST

176:    test:
177:       suffix: 1
178:       requires: !single

180:    test:
181:       suffix: 2
182:       args: -rank 40
183:       requires: !single

185: TEST*/