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