Actual source code: ex32.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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 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;
 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:   PetscOptionsGetInt(NULL,NULL,"-rank",&rank,NULL);
 38:   PetscPrintf(PETSC_COMM_WORLD,"\nLyapunov equation, N=%D (%Dx%D grid)\n\n",N,n,m);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

156:   LMEGetErrorEstimate(lme,&errest);
157:   PetscPrintf(PETSC_COMM_WORLD," Error estimate reported by the solver: %.4g\n",(double)errest);
158:   LMEComputeError(lme,&error);
159:   PetscPrintf(PETSC_COMM_WORLD," Computed residual norm: %.4g\n\n",(double)error);

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

171: /*TEST

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

177: TEST*/