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