1/*******************************************************************************
2 * Copyright (c) 2018, College of William & Mary
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *     * Redistributions of source code must retain the above copyright
8 *       notice, this list of conditions and the following disclaimer.
9 *     * Redistributions in binary form must reproduce the above copyright
10 *       notice, this list of conditions and the following disclaimer in the
11 *       documentation and/or other materials provided with the distribution.
12 *     * Neither the name of the College of William & Mary nor the
13 *       names of its contributors may be used to endorse or promote products
14 *       derived from this software without specific prior written permission.
15 *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19 * DISCLAIMED. IN NO EVENT SHALL THE COLLEGE OF WILLIAM & MARY BE LIABLE FOR ANY
20 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 *
27 * PRIMME: https://github.com/primme/primme
28 * Contact: Andreas Stathopoulos, a n d r e a s _at_ c s . w m . e d u
29 *******************************************************************************
30 *
31 *  Example to compute the k largest eigenvalues in a 1-D Laplacian matrix.
32 *
33 ******************************************************************************/
34
35#include <stdlib.h>
36#include <stdio.h>
37#include <math.h>
38ifdef(`USE_COMPLEX', ifdef(`USE_COMPLEX_CXX', ``#include <complex>'', ``#include <complex.h>''))
39ifdef(`USE_PETSC', ``#include <petscpc.h>
40#include <petscmat.h>
41'')dnl
42#include "primme.h"   /* header file is required to run primme */
43define(`PRIMME_REAL', ifdef(`USE_PETSC', `PetscReal', ifdef(`USE_FLOAT', `float', `double')))dnl
44define(`PRIMME_SCALAR', ifdef(`USE_PETSC', `PetscScalar', ifdef(`USE_COMPLEX', ifdef(`USE_COMPLEX_CXX', `std::complex<PRIMME_REAL>', `complex PRIMME_REAL'), `PRIMME_REAL')))dnl
45ifdef(`USE_PETSC', `
46PetscErrorCode generateLaplacian1D(int n, Mat *A);
47void PETScMatvec(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *ierr);
48void ApplyPCPrecPETSC(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *ierr);
49static void par_GlobalSum(void *sendBuf, void *recvBuf, int *count,
50                         primme_params *primme, int *ierr);
51', `
52void LaplacianMatrixMatvec(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *ierr);
53void LaplacianApplyPreconditioner(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *ierr);
54')dnl
55
56int main (int argc, char *argv[]) {
57
58   /* Solver arrays and parameters */
59   PRIMME_REAL *evals;    /* Array with the computed eigenvalues */
60   PRIMME_REAL *rnorms;   /* Array with the computed eigenpairs residual norms */
61   PRIMME_SCALAR *evecs;    /* Array with the computed eigenvectors;
62                        first vector starts in evecs[0],
63                        second vector starts in evecs[primme.n],
64                        third vector starts in evecs[primme.n*2]...  */
65   primme_params primme;
66                     /* PRIMME configuration struct */
67ifdef(`ADVANCED', `   double targetShifts[1];
68')dnl
69
70   /* Other miscellaneous items */
71   int ret;
72   int i;
73ifdef(`USE_PETSC', `   Mat A; /* problem matrix */
74   PC pc;            /* preconditioner */
75   PetscErrorCode ierr;
76   PetscInt n, nLocal;
77   MPI_Comm comm;
78
79   PetscInitialize(&argc, &argv, NULL, NULL);
80
81')dnl
82
83   /* Set default values in PRIMME configuration struct */
84   primme_initialize(&primme);
85
86   /* Set problem matrix */
87ifdef(`USE_PETSC', `   ierr = generateLaplacian1D(100, &A); CHKERRQ(ierr);
88   primme.matrix = &A;
89   primme.matrixMatvec = PETScMatvec;
90', `   primme.matrixMatvec = LaplacianMatrixMatvec;
91')dnl
92                           /* Function that implements the matrix-vector product
93                              A*x for solving the problem A*x = l*x */
94
95   /* Set problem parameters */
96ifdef(`USE_PETSC', `   ierr = MatGetSize(A, &n, NULL); CHKERRQ(ierr);
97   primme.n = (PRIMME_INT)n;',
98`   primme.n = 100;') /* set problem dimension */
99   primme.numEvals = 10;   /* Number of wanted eigenpairs */
100   primme.eps = ifdef(`USE_PETSC', 1e-5, 1e-9);      /* ||r|| <= eps * ||matrix|| */
101   primme.target = primme_smallest;
102                           /* Wanted the smallest eigenvalues */
103
104   /* Set preconditioner (optional) */
105ifdef(`USE_PETSC', `   ierr = PCCreate(PETSC_COMM_WORLD, &pc); CHKERRQ(ierr);
106   ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr);
107   ierr = PCSetOperators(pc, A, A); CHKERRQ(ierr);
108   ierr = PCSetFromOptions(pc); CHKERRQ(ierr);
109   ierr = PCSetUp(pc); CHKERRQ(ierr);
110   primme.preconditioner = &pc;
111   primme.applyPreconditioner = ApplyPCPrecPETSC;
112', `   primme.applyPreconditioner = LaplacianApplyPreconditioner;
113')dnl
114   primme.correctionParams.precondition = 1;
115
116   /* Set advanced parameters if you know what are you doing (optional) */
117   /*
118   primme.maxBasisSize = 14;
119   primme.minRestartSize = 4;
120   primme.maxBlockSize = 1;
121   primme.maxMatvecs = 1000;
122   */
123
124   /* Set method to solve the problem */
125   primme_set_method(PRIMME_DYNAMIC, &primme);
126   /* DYNAMIC uses a runtime heuristic to choose the fastest method between
127       PRIMME_DEFAULT_MIN_TIME and PRIMME_DEFAULT_MIN_MATVECS. But you can
128       set another method, such as PRIMME_LOBPCG_OrthoBasis_Window, directly */
129
130ifdef(`USE_PETSC', `   /* Set parallel parameters */
131   ierr = MatGetLocalSize(A, &nLocal, NULL); CHKERRQ(ierr);
132   primme.nLocal = (PRIMME_INT)nLocal;
133   comm = PETSC_COMM_WORLD;
134   primme.commInfo = &comm;
135   MPI_Comm_size(comm, &primme.numProcs);
136   MPI_Comm_rank(comm, &primme.procID);
137   primme.globalSumReal = par_GlobalSum;
138
139')dnl
140   /* Display PRIMME configuration struct (optional) */
141ifdef(`USE_PETSC', `   if (primme.procID == 0) /* Reports process with ID 0 */
142   ')   primme_display_params(primme);
143
144   /* Allocate space for converged Ritz values and residual norms */ifdef(`USE_COMPLEX_CXX', `
145   evals = new PRIMME_REAL[primme.numEvals];
146   evecs = new PRIMME_SCALAR[primme.n*primme.numEvals];
147   rnorms = new PRIMME_REAL[primme.numEvals];',`
148   evals = (PRIMME_REAL*)malloc(primme.numEvals*sizeof(PRIMME_REAL));
149   evecs = (PRIMME_SCALAR*)malloc(primme.n*primme.numEvals*sizeof(PRIMME_SCALAR));
150   rnorms = (PRIMME_REAL*)malloc(primme.numEvals*sizeof(PRIMME_REAL));')
151
152define(`CALL_PRIMME', `   /* Call primme  */
153ifdef(`USE_PETSC', ``#if defined(PETSC_USE_COMPLEX) && defined(PETSC_USE_REAL_SINGLE)
154   ret = cprimme(evals, evecs, rnorms, &primme);
155#elif defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE)
156   ret = zprimme(evals, evecs, rnorms, &primme);
157#elif !defined(PETSC_USE_COMPLEX) && defined(PETSC_USE_REAL_SINGLE)
158   ret = sprimme(evals, evecs, rnorms, &primme);
159#elif !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE)
160   ret = dprimme(evals, evecs, rnorms, &primme);
161#endif
162'',
163`   ret = ifdef(`USE_COMPLEX', ifdef(`USE_FLOAT',`c',`z'), ifdef(`USE_FLOAT',`s',`d'))primme(evals, evecs, rnorms, &primme);
164')dnl
165
166   if (ret != 0) {
167      fprintf(primme.outputFile,
168         "Error: primme returned with nonzero exit status: %d \n",ret);
169      return -1;
170   }
171
172ifdef(`USE_PETSC', ``   if (primme.procID == 0) { /* Reports process with ID 0 */
173' define(sp, `   ')', `define(sp, `')')dnl
174   sp()/* Reporting (optional) */
175   sp()for (i=0; i < primme.initSize; i++) {
176   sp()   fprintf(primme.outputFile, "Eval[%d]: %-22.15E rnorm: %-22.15E\n", i+1,
177   sp()      evals[i], rnorms[i]);
178   sp()}
179   sp()fprintf(primme.outputFile, " %d eigenpairs converged\n", primme.initSize);
180   sp()fprintf(primme.outputFile, "Tolerance : %-22.15E\n",
181   sp()                                                      primme.aNorm*primme.eps);
182   sp()fprintf(primme.outputFile, "Iterations: %-" PRIMME_INT_P "\n",
183   sp()                                              primme.stats.numOuterIterations);
184   sp()fprintf(primme.outputFile, "Restarts  : %-" PRIMME_INT_P "\n", primme.stats.numRestarts);
185   sp()fprintf(primme.outputFile, "Matvecs   : %-" PRIMME_INT_P "\n", primme.stats.numMatvecs);
186   sp()fprintf(primme.outputFile, "Preconds  : %-" PRIMME_INT_P "\n", primme.stats.numPreconds);
187   sp()if (primme.stats.lockingIssue) {
188   sp()   fprintf(primme.outputFile, "\nA locking problem has occurred.\n");
189   sp()   fprintf(primme.outputFile,
190   sp()      "Some eigenpairs do not have a residual norm less than the tolerance.\n");
191   sp()   fprintf(primme.outputFile,
192   sp()      "However, the subspace of evecs is accurate to the required tolerance.\n");
193   sp()}
194
195   sp()switch (primme.dynamicMethodSwitch) {
196   sp()   case -1: fprintf(primme.outputFile,
197   sp()         "Recommended method for next run: DEFAULT_MIN_MATVECS\n"); break;
198   sp()   case -2: fprintf(primme.outputFile,
199   sp()         "Recommended method for next run: DEFAULT_MIN_TIME\n"); break;
200   sp()   case -3: fprintf(primme.outputFile,
201   sp()         "Recommended method for next run: DYNAMIC (close call)\n"); break;
202   sp()}
203ifdef(`USE_PETSC', `   }
204')dnl
205')dnl end of CALL_PRIMME
206CALL_PRIMME
207ifdef(`ADVANCED', `
208   /* Note that d/zprimme can be called more than once before call primme_free. */
209   /* Find the 5 eigenpairs closest to .5 */
210   primme.numTargetShifts = 1;
211   targetShifts[0] = .5;
212   primme.targetShifts = targetShifts;
213   primme.target = primme_closest_abs;
214   primme.numEvals = 5;
215   primme.initSize = 0; /* primme.initSize may be not zero after a d/zprimme;
216                           so set it to zero to avoid the already converged eigenvectors
217                           being used as initial vectors. */
218
219CALL_PRIMME
220
221   /* Perturb the 5 approximate eigenvectors in evecs and used them as initial solution.
222      This time the solver should converge faster than the last one. */
223   for (i=0; i<primme.n*5; i++)
224      evecs[i] += rand()/(PRIMME_REAL)RAND_MAX*1e-4;
225   primme.initSize = 5;
226   primme.numEvals = 5;
227
228CALL_PRIMME
229
230   /* Find the next 5 eigenpairs closest to .5 */
231   primme.initSize = 0;
232   primme.numEvals = 5;
233   primme.numOrthoConst = 5; /* solver will find solutions orthogonal to the already
234                                5 approximate eigenvectors in evecs */
235
236CALL_PRIMME
237')dnl
238   primme_free(&primme);ifdef(`USE_COMPLEX_CXX', `
239   delete [] evals;
240   delete [] evecs;
241   delete [] rnorms;',`
242   free(evals);
243   free(evecs);
244   free(rnorms);')
245
246ifdef(`USE_PETSC', `   ierr = PetscFinalize(); CHKERRQ(ierr);
247
248')dnl
249  return(0);
250}
251
252/* 1-D Laplacian block matrix-vector product, Y = A * X, where
253
254   - X, input dense matrix of size primme.n x blockSize;
255   - Y, output dense matrix of size primme.n x blockSize;
256   - A, tridiagonal square matrix of dimension primme.n with this form:
257
258        [ 2 -1  0  0  0 ... ]
259        [-1  2 -1  0  0 ... ]
260        [ 0 -1  2 -1  0 ... ]
261         ...
262*/
263ifdef(`USE_PETSC', `
264PetscErrorCode generateLaplacian1D(int n, Mat *A) {
265   PetscScalar    value[3] = {-1.0, 2.0, -1.0};
266   PetscInt       i,Istart,Iend,col[3];
267   PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
268   PetscErrorCode ierr;
269
270   PetscFunctionBegin;
271
272   ierr = MatCreate(PETSC_COMM_WORLD, A); CHKERRQ(ierr);
273   ierr = MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, n); CHKERRQ(ierr);
274   ierr = MatSetFromOptions(*A); CHKERRQ(ierr);
275   ierr = MatSetUp(*A); CHKERRQ(ierr);
276
277   ierr = MatGetOwnershipRange(*A, &Istart, &Iend); CHKERRQ(ierr);
278   if (Istart == 0) FirstBlock = PETSC_TRUE;
279   if (Iend == n) LastBlock = PETSC_TRUE;
280   for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
281      col[0]=i-1; col[1]=i; col[2]=i+1;
282      ierr = MatSetValues(*A, 1, &i, 3, col, value, INSERT_VALUES); CHKERRQ(ierr);
283   }
284   if (LastBlock) {
285      i=n-1; col[0]=n-2; col[1]=n-1;
286      ierr = MatSetValues(*A, 1, &i, 2, col, value, INSERT_VALUES); CHKERRQ(ierr);
287   }
288   if (FirstBlock) {
289      i=0; col[0]=0; col[1]=1; value[0]=2.0; value[1]=-1.0;
290      ierr = MatSetValues(*A, 1, &i, 2, col, value, INSERT_VALUES); CHKERRQ(ierr);
291   }
292
293   ierr = MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
294   ierr = MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
295
296   PetscFunctionReturn(0);
297}
298
299void PETScMatvec(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *err) {
300   int i;
301   Mat *matrix;
302   Vec xvec, yvec;
303   PetscErrorCode ierr;
304
305   matrix = (Mat *)primme->matrix;
306
307   ierr = MatCreateVecs(*matrix, &xvec, &yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
308   for (i=0; i<*blockSize; i++) {
309      ierr = VecPlaceArray(xvec, ((PetscScalar*)x) + *ldx*i); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
310      ierr = VecPlaceArray(yvec, ((PetscScalar*)y) + *ldy*i); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
311      ierr = MatMult(*matrix, xvec, yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
312      ierr = VecResetArray(xvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
313      ierr = VecResetArray(yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
314   }
315   ierr = VecDestroy(&xvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
316   ierr = VecDestroy(&yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
317   *err = 0;
318}
319', `
320void LaplacianMatrixMatvec(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *err) {
321
322   int i;            /* vector index, from 0 to *blockSize-1*/
323   int row;          /* Laplacian matrix row index, from 0 to matrix dimension */
324   PRIMME_SCALAR *xvec;     /* pointer to i-th input vector x */
325   PRIMME_SCALAR *yvec;     /* pointer to i-th output vector y */
326
327   for (i=0; i<*blockSize; i++) {
328      xvec = (PRIMME_SCALAR *)x + *ldx*i;
329      yvec = (PRIMME_SCALAR *)y + *ldy*i;
330      for (row=0; row<primme->n; row++) {
331         yvec[row] = 0.0;
332         if (row-1 >= 0) yvec[row] += -1.0*xvec[row-1];
333         yvec[row] += 2.0*xvec[row];
334         if (row+1 < primme->n) yvec[row] += -1.0*xvec[row+1];
335      }
336   }
337   *err = 0;
338}
339')dnl
340
341/* This performs Y = M^{-1} * X, where
342
343   - X, input dense matrix of size primme.n x blockSize;
344   - Y, output dense matrix of size primme.n x blockSize;
345   - M, diagonal square matrix of dimension primme.n with 2 in the diagonal.
346*/
347ifdef(`USE_PETSC', `
348void ApplyPCPrecPETSC(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *err) {
349   int i;
350   Mat *matrix;
351   PC *pc;
352   Vec xvec, yvec;
353   PetscErrorCode ierr;
354
355   matrix = (Mat *)primme->matrix;
356   pc = (PC *)primme->preconditioner;
357
358   ierr = MatCreateVecs(*matrix, &xvec, &yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
359   for (i=0; i<*blockSize; i++) {
360      ierr = VecPlaceArray(xvec, ((PetscScalar*)x) + *ldx*i); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
361      ierr = VecPlaceArray(yvec, ((PetscScalar*)y) + *ldy*i); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
362      ierr = PCApply(*pc, xvec, yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
363      ierr = VecResetArray(xvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
364      ierr = VecResetArray(yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
365   }
366   ierr = VecDestroy(&xvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
367   ierr = VecDestroy(&yvec); CHKERRABORT(*(MPI_Comm*)primme->commInfo, ierr);
368}
369
370static void par_GlobalSum(void *sendBuf, void *recvBuf, int *count,
371                         primme_params *primme, int *ierr) {
372   MPI_Comm communicator = *(MPI_Comm *) primme->commInfo;
373
374   if (sendBuf == recvBuf) {
375     *ierr = MPI_Allreduce(MPI_IN_PLACE, recvBuf, *count, MPIU_REAL, MPI_SUM, communicator) != MPI_SUCCESS;
376   } else {
377     *ierr = MPI_Allreduce(sendBuf, recvBuf, *count, MPIU_REAL, MPI_SUM, communicator) != MPI_SUCCESS;
378   }
379}
380', `
381void LaplacianApplyPreconditioner(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *ierr) {
382
383   int i;            /* vector index, from 0 to *blockSize-1*/
384   int row;          /* Laplacian matrix row index, from 0 to matrix dimension */
385   PRIMME_SCALAR *xvec;     /* pointer to i-th input vector x */
386   PRIMME_SCALAR *yvec;     /* pointer to i-th output vector y */
387
388   for (i=0; i<*blockSize; i++) {
389      xvec = (PRIMME_SCALAR *)x + *ldx*i;
390      yvec = (PRIMME_SCALAR *)y + *ldy*i;
391      for (row=0; row<primme->n; row++) {
392         yvec[row] = xvec[row]/2.;
393      }
394   }
395   *ierr = 0;
396}
397')dnl
398