1C******************************************************************************* 2C Copyright (c) 2018, College of William & Mary 3C All rights reserved. 4C 5C Redistribution and use in source and binary forms, with or without 6C modification, are permitted provided that the following conditions are met: 7C * Redistributions of source code must retain the above copyright 8C notice, this list of conditions and the following disclaimer. 9C * Redistributions in binary form must reproduce the above copyright 10C notice, this list of conditions and the following disclaimer in the 11C documentation and/or other materials provided with the distribution. 12C * Neither the name of the College of William & Mary nor the 13C names of its contributors may be used to endorse or promote products 14C derived from this software without specific prior written permission. 15C 16C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 17C ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 18C WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 19C DISCLAIMED. IN NO EVENT SHALL THE COLLEGE OF WILLIAM & MARY BE LIABLE FOR ANY 20C DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 21C (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 22C LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 23C ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 24C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 25C SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26C 27C PRIMME: https://github.com/primme/primme 28C 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 Program primmeF77Example 36!----------------------------------------------------------------------- 37 implicit none 38 include 'primme_f77.h' 39! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 40! Pointer to the PRIMME data structure used internally by PRIMME 41! 42! Note that for 64 bit systems, pointers are 8 bytes so use: 43 integer*8 primme 44! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 45! Problem setup 46! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47 48 ! Solver Parameters 49 integer*8 n,NUMEmax,BASISmax,BLOCKmax,maxMatvecs, 50 : printLevel, method, whichEvals, numTargetShifts 51 real*8 ETOL 52 53 parameter ( 54 : n = 100, 55 : BASISmax = 12, 56 : NUMEmax = 5, 57 : BLOCKmax = 1, 58 : maxMatvecs = 300000, 59 : ETOL = 1.0D-12, 60 : printLevel = 5, 61 : whichEvals = primme_smallest, 62 : numTargetShifts = 2, 63 : method = PRIMME_DYNAMIC 64 : ) 65 real*8 TargetShifts(numTargetShifts) 66 data TargetShifts /3.0D0, 5.1D0/ 67 68 external MV, ApplyPrecon 69 70! Eigenvalues, eigenvectors, and their residual norms 71! 72 real*8 evals(NUMEmax) 73 real*8 rnorms(NUMEmax) 74 real*8 evecs(n*NUMEmax) 75 76! Other vars 77! 78 integer i,ierr 79 real*8 epsil, aNorm 80 integer*8 numIts, numMatvecs 81 82!----------------------------------------------------------------------- 83! Start executable 84!----------------------------------------------------------------------- 85! 86! ---------------------------------------------------------------- 87! Initialize PRIMME 88! ---------------------------------------------------------------- 89! 90 call primme_initialize_f77(primme) 91 92! Set a few basic solver parameters 93 call primme_set_member_f77(primme, PRIMME_n, n, ierr) 94 call primme_set_member_f77(primme, PRIMME_numEvals, NUMEmax, 95 : ierr) 96 call primme_set_member_f77(primme, PRIMME_eps, ETOL, ierr) 97 call primme_set_member_f77(primme, PRIMME_target, 98 : whichEvals, ierr) 99 call primme_set_member_f77(primme, PRIMME_numTargetShifts, 100 : numTargetShifts, ierr) 101 call primme_set_member_f77(primme, PRIMME_targetShifts, 102 : TargetShifts, ierr) 103 104! Set matvec 105 call primme_set_member_f77(primme, PRIMME_matrixMatvec, 106 : MV, ierr) 107 108 109! Set preconditioner (optional) 110 call primme_set_member_f77(primme, 111 : PRIMME_applyPreconditioner, ApplyPrecon, ierr) 112 call primme_set_member_f77(primme, 113 : PRIMME_correctionParams_precondition, 0, ierr) 114! 115! Set a few other solver parameters (optional) 116! 117 call primme_set_member_f77(primme, PRIMME_maxBasisSize, 118 : BASISmax, ierr) 119 call primme_set_member_f77(primme, PRIMME_maxBlockSize, 120 : BLOCKmax, ierr) 121 call primme_set_member_f77(primme, PRIMME_printLevel, 122 : printLevel, ierr) 123 call primme_set_member_f77(primme, PRIMME_maxMatvecs, 124 : maxMatvecs, ierr) 125! 126! Set the method to be used (after n, numEvals, and precondition have 127! been set. Also after basisSize is set, if desired.) 128 call primme_set_method_f77(primme, method, ierr) 129 130 if (ierr .lt. 0) 131 : write(*,*) 'No preset method. Using custom settings' 132 133! ---------------------------------------------------------------- 134! Display what parameters are used 135! ---------------------------------------------------------------- 136 137 call primme_display_params_f77(primme) 138 139! ---------------------------------------------------------------- 140! Calling the PRIMME solver 141! ---------------------------------------------------------------- 142 143 call dprimme_f77(evals, evecs, rnorms, primme, ierr) 144 if (ierr.ne.0) then 145 stop 1 146 endif 147 148! ---------------------------------------------------------------- 149! Reporting results 150 151 if (ierr.eq.0) then 152 print *, 'PRIMME has returned successfully' 153 else 154 print *, 'PRIMME returned with error: ', ierr 155 endif 156 157! 158! Example of obtaining primme members from the driver: 159! NOTE: don't use primme_get_member_f77, which can only be used in a callback 160! 161 call primmetop_get_member_f77(primme, PRIMME_eps, epsil, 162 : ierr) 163 call primmetop_get_member_f77(primme, PRIMME_aNorm, 164 : aNorm, ierr) 165 call primmetop_get_member_f77(primme, 166 : PRIMME_stats_numOuterIterations, numIts, ierr) 167 call primmetop_get_member_f77(primme, 168 : PRIMME_stats_numMatvecs, numMatvecs, ierr) 169 print '(A,E8.2,/,A,e12.5,/,A,I8,/,A,I8)', 170 : 'Tolerance used: ',epsil, 171 : 'Estimated norm(A):',aNorm, 172 : 'Iterations: ',numIts, 173 : 'Matvecs: ',numMatvecs 174! 175! Reporting of evals and residuals 176! 177 do i = 1, numemax 178 write (*, 9000) i, evals(i),rnorms(i) 179 enddo 180 9000 FORMAT (1x,'E(',i1,') = ',G24.16,4x, 181 & 'residual norm =', E12.4) 182 call primme_free_f77(primme) 183 stop 184 write(0,*) 'ERROR! No data in the file' 185 stop 186 end 187!----------------------------------------------------------------------- 188! Supporting subroutines 189!----------------------------------------------------------------------- 190! ---------------------------------------------------------------- 191 192 193! 1-D Laplacian block matrix-vector product, Y = A * X, where 194! 195! - X, input dense matrix of size primme.n x blockSize; 196! - Y, output dense matrix of size primme.n x blockSize; 197! - A, tridiagonal square matrix of dimension primme.n with this form: 198! 199! 2 -1 0 0 0 ... 200! -1 2 -1 0 0 ... 201! 0 -1 2 -1 0 ... 202! ... 203! 204 subroutine MV(x,ldx,y,ldy,k,primme,ierr) 205! ---------------------------------------------------------------- 206 implicit none 207 include 'primme_f77.h' 208 integer*8 ldx, ldy 209 real*8 x(ldx,*), y(ldy,*) 210 integer*8 primme 211 integer*8 n, i 212 integer k,ierr,j 213 call primme_get_member_f77(primme, PRIMME_n, n, ierr) 214 do j=1,k 215 do i=1,n 216 y(i,j) = 0 217 if (i.ge.2) then 218 y(i,j) = y(i,j) - x(i-1,j) 219 endif 220 y(i,j) = y(i,j) + 2.*x(i,j) 221 if (i.le.n-1) then 222 y(i,j) = y(i,j) - x(i+1,j) 223 endif 224 enddo 225 enddo 226 ierr = 0 227 end 228 229! This performs Y = M^{-1} * X, where 230! 231! - X, input dense matrix of size primme.n x blockSize; 232! - Y, output dense matrix of size primme.n x blockSize; 233! - M, diagonal square matrix of dimension primme.n with 2 in the diagonal. 234! 235 subroutine ApplyPrecon(x,ldx,y,ldy,k,primme, ierr) 236! ---------------------------------------------------------------- 237 implicit none 238 include 'primme_f77.h' 239 integer*8 ldx, ldy 240 real*8 x(ldx,*), y(ldy,*) 241 integer*8 primme 242 integer*8 n, i 243 integer k,ierr,j 244 call primme_get_member_f77(primme, PRIMME_n, n, ierr) 245 do j=1,k 246 do i=1,n 247 y(i,j) = x(i,j)/2.0 248 enddo 249 enddo 250 ierr = 0 251 end 252