Actual source code: ex32.c

slepc-main 2024-11-15
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,NULL,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(MatGetOwnershipRange(A,&Istart,&Iend));
 49:   for (II=Istart;II<Iend;II++) {
 50:     i = II/n; j = II-i*n;
 51:     if (i>0) PetscCall(MatSetValue(A,II,II-n,1.0,INSERT_VALUES));
 52:     if (i<m-1) PetscCall(MatSetValue(A,II,II+n,1.0,INSERT_VALUES));
 53:     if (j>0) PetscCall(MatSetValue(A,II,II-1,1.0,INSERT_VALUES));
 54:     if (j<n-1) PetscCall(MatSetValue(A,II,II+1,1.0,INSERT_VALUES));
 55:     PetscCall(MatSetValue(A,II,II,-4.0-sigma,INSERT_VALUES));
 56:   }
 57:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
 58:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));

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

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

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

 89:   /*
 90:      Set the type of equation
 91:   */
 92:   PetscCall(LMESetProblemType(lme,LME_LYAPUNOV));

 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:   PetscCall(LMESetCoefficients(lme,A,NULL,NULL,NULL));
100:   PetscCall(LMESetRHS(lme,C));

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

115:   /*
116:      (Optional) Set other solver options
117:   */
118:   PetscCall(LMESetTolerances(lme,1e-07,PETSC_CURRENT));

120:   /*
121:      Set solver parameters at runtime
122:   */
123:   PetscCall(LMESetFromOptions(lme));

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:                    Solve the matrix equation, A*X+X*A'=-C
127:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

129:   PetscCall(LMESolve(lme));
130:   PetscCall(LMEGetConvergedReason(lme,&reason));
131:   PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");

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

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

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:                         Compute residual error
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

171: /*TEST

173:    test:
174:       suffix: 1
175:       requires: !single

177:    test:
178:       suffix: 2
179:       args: -rank 40
180:       requires: !single

182: TEST*/