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 largest eigenvalues in a diagonal matrix in 32 * parallel using MPI. 33 * 34 ******************************************************************************/ 35 36 #include <stdlib.h> 37 #include <stdio.h> 38 #include <math.h> 39 #include <mpi.h> 40 #include <assert.h> 41 42 #include "primme.h" /* header file is required to run primme */ 43 44 void DiagonalMatrixMatvec(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *ierr); 45 static void par_GlobalSum(void *sendBuf, void *recvBuf, int *count, 46 primme_params *primme, int *ierr); 47 48 #ifndef min 49 # define min(a, b) ((a) < (b) ? (a) : (b)) 50 #endif 51 52 int main (int argc, char *argv[]) { 53 54 /* Solver arrays and parameters */ 55 float *evals; /* Array with the computed eigenvalues */ 56 float *rnorms; /* Array with the computed eigenpairs residual norms */ 57 float *evecs; /* Array with the computed eigenvectors; 58 first vector starts in evecs[0], 59 second vector starts in evecs[primme.n], 60 third vector starts in evecs[primme.n*2]... */ 61 primme_params primme; 62 /* PRIMME configuration struct */ 63 64 /* Other miscellaneous items */ 65 int ret; 66 int i; 67 68 /* Initialize the infrastructure necessary for communication */ 69 MPI_Init(&argc, &argv); 70 71 /* Set default values in PRIMME configuration struct */ 72 primme_initialize(&primme); 73 74 /* Set problem matrix */ 75 primme.matrixMatvec = DiagonalMatrixMatvec; 76 /* Function that implements the matrix-vector product 77 A*x for solving the problem A*x = l*x */ 78 79 /* Set problem parameters */ 80 primme.n = 1000; /* set problem dimension */ 81 primme.numEvals = 1000; /* Number of wanted eigenpairs */ 82 primme.eps = .1; /* ||r|| <= eps * ||matrix|| */ 83 primme.target = primme_largest; 84 /* Wanted the smallest eigenvalues */ 85 86 /* Set advanced parameters if you know what are you doing (optional) */ 87 /* 88 primme.maxBasisSize = 14; 89 primme.minRestartSize = 4; 90 primme.maxBlockSize = 1; 91 primme.maxMatvecs = 1000; 92 */ 93 94 /* Set method to solve the problem */ 95 primme_set_method(PRIMME_DEFAULT_MIN_MATVECS, &primme); 96 /* DYNAMIC uses a runtime heuristic to choose the fastest method between 97 PRIMME_DEFAULT_MIN_TIME and PRIMME_DEFAULT_MIN_MATVECS. But you can 98 set another method, such as PRIMME_LOBPCG_OrthoBasis_Window, directly */ 99 100 /* Set parallel parameters */ 101 MPI_Comm comm = MPI_COMM_WORLD; 102 MPI_Comm_size(comm, &primme.numProcs); 103 MPI_Comm_rank(comm, &primme.procID); 104 primme.commInfo = &comm; /* User-defined member to pass the communicator to 105 globalSumReal and broadcastReal */ 106 /* In this example, the matrix is distributed by rows, and the first 107 * processes may have an extra row in order to distribute the remaining rows 108 * n % numProcs */ 109 PRIMME_INT nLocal = primme.n / primme.numProcs + 110 (primme.n % primme.numProcs > primme.procID ? 1 : 0); 111 primme.nLocal = nLocal; /* Number of local rows */ 112 primme.globalSumReal = par_GlobalSum; 113 114 /* Display PRIMME configuration struct (optional) */ 115 if (primme.procID == 0) primme_display_params(primme); 116 117 /* Allocate space for converged Ritz values and residual norms */ 118 evals = (float*)malloc(primme.numEvals*sizeof(float)); 119 evecs = (float*)malloc(primme.n*primme.numEvals*sizeof(float)); 120 rnorms = (float*)malloc(primme.numEvals*sizeof(float)); 121 122 /* Call primme */ 123 ret = sprimme(evals, evecs, rnorms, &primme); 124 125 if (primme.procID == 0) { 126 if (ret != 0) { 127 fprintf(primme.outputFile, 128 "Error: primme returned with nonzero exit status: %d \n",ret); 129 return -1; 130 } 131 132 /* Reporting (optional) */ 133 for (i=0; i < primme.initSize; i++) { 134 fprintf(primme.outputFile, "Eval[%d]: %-22.15E rnorm: %-22.15E\n", i+1, 135 evals[i], rnorms[i]); 136 } 137 fprintf(primme.outputFile, " %d eigenpairs converged\n", primme.initSize); 138 fprintf(primme.outputFile, "Tolerance : %-22.15E\n", 139 primme.aNorm*primme.eps); 140 fprintf(primme.outputFile, "Iterations: %-" PRIMME_INT_P "\n", 141 primme.stats.numOuterIterations); 142 fprintf(primme.outputFile, "Restarts : %-" PRIMME_INT_P "\n", primme.stats.numRestarts); 143 fprintf(primme.outputFile, "Matvecs : %-" PRIMME_INT_P "\n", primme.stats.numMatvecs); 144 fprintf(primme.outputFile, "Preconds : %-" PRIMME_INT_P "\n", primme.stats.numPreconds); 145 fprintf(primme.outputFile, "Orthogonalization Time : %g\n", primme.stats.timeOrtho); 146 fprintf(primme.outputFile, "Matvec Time : %g\n", primme.stats.timeMatvec); 147 fprintf(primme.outputFile, "GlobalSum Time : %g\n", primme.stats.timeGlobalSum); 148 fprintf(primme.outputFile, "Broadcast Time : %g\n", primme.stats.timeBroadcast); 149 fprintf(primme.outputFile, "Total Time : %g\n", primme.stats.elapsedTime); 150 if (primme.stats.lockingIssue) { 151 fprintf(primme.outputFile, "\nA locking problem has occurred.\n"); 152 fprintf(primme.outputFile, 153 "Some eigenpairs do not have a residual norm less than the tolerance.\n"); 154 fprintf(primme.outputFile, 155 "However, the subspace of evecs is accurate to the required tolerance.\n"); 156 } 157 158 switch (primme.dynamicMethodSwitch) { 159 case -1: fprintf(primme.outputFile, 160 "Recommended method for next run: DEFAULT_MIN_MATVECS\n"); break; 161 case -2: fprintf(primme.outputFile, 162 "Recommended method for next run: DEFAULT_MIN_TIME\n"); break; 163 case -3: fprintf(primme.outputFile, 164 "Recommended method for next run: DYNAMIC (close call)\n"); break; 165 } 166 } 167 168 primme_free(&primme); 169 free(evals); 170 free(evecs); 171 free(rnorms); 172 173 /* Tear down the communication infrastructure */ 174 MPI_Finalize(); 175 176 return(0); 177 } 178 179 /* Diagonal block matrix-vector product, Y = A * X, where 180 181 - X, input dense matrix of size primme.n x blockSize; 182 - Y, output dense matrix of size primme.n x blockSize; 183 - A, tridiagonal square matrix of dimension primme.n with this form: 184 185 */ 186 187 void DiagonalMatrixMatvec(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize, primme_params *primme, int *err) { 188 189 int i; /* vector index, from 0 to *blockSize-1*/ 190 int row; /* local matrix row index, from 0 to nLocal */ 191 /* In this example, row0 is the global index of the first local row */ 192 int row0 = primme->n / primme->numProcs * primme->procID + 193 min(primme->n % primme->numProcs, primme->procID); 194 float *xvec; /* pointer to i-th input vector x */ 195 float *yvec; /* pointer to i-th output vector y */ 196 197 for (i=0; i<*blockSize; i++) { 198 xvec = (float *)x + *ldx*i; 199 yvec = (float *)y + *ldy*i; 200 for (row = 0; row < primme->nLocal; row++) { 201 /* The diagonal matrix has the spectrum of a Laplacial */ 202 float v = sin(M_PI * (row + row0 + 1) / 2.0 / (primme->n + 1)); 203 yvec[row] = 4. * v * v * xvec[row]; 204 } 205 } 206 *err = 0; 207 } 208 209 static void par_GlobalSum(void *sendBuf, void *recvBuf, int *count, 210 primme_params *primme, int *ierr) { 211 MPI_Comm communicator = *(MPI_Comm *) primme->commInfo; 212 213 if (sendBuf == recvBuf) { 214 *ierr = MPI_Allreduce(MPI_IN_PLACE, recvBuf, *count, MPI_FLOAT, MPI_SUM, communicator) != MPI_SUCCESS; 215 } else { 216 *ierr = MPI_Allreduce(sendBuf, recvBuf, *count, MPI_FLOAT, MPI_SUM, communicator) != MPI_SUCCESS; 217 } 218 } 219