source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/examples/dsdrv6.f @ 3225

Last change on this file since 3225 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 17.5 KB
Line 
1      program dsdrv6 
2c
3c     Program to illustrate the idea of reverse communication
4c     in Cayley mode for a generalized symmetric eigenvalue 
5c     problem.  The following program uses the two LAPACK subroutines 
6c     dgttrf.f and dgttrs.f to factor and solve a tridiagonal system of 
7c     equations.
8c
9c     We implement example six of ex-sym.doc in DOCUMENTS directory
10c
11c\Example-6
12c     ... Suppose we want to solve A*x = lambda*M*x in inverse mode,
13c         where A and M are obtained by the finite element of the
14c         1-dimensional discrete Laplacian
15c                             d^2u / dx^2
16c         on the interval [0,1] with zero Dirichlet boundary condition
17c         using piecewise linear elements.
18c
19c     ... OP = (inv[A-sigma*M])*(A+sigma*M)  and  B = M.
20c
21c     ... Use mode 5 of DSAUPD.
22c
23c\BeginLib
24c
25c\References:
26c  1. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
27c     Algorithm for Solving Sparse Symmetric Generalized Eigenproblems", 
28c     SIAM J. Matr. Anal. Apps.,  January (1993).
29c
30c\Routines called:
31c     dsaupd  ARPACK reverse communication interface routine.
32c     dseupd  ARPACK routine that returns Ritz values and (optionally)
33c             Ritz vectors.
34c     dgttrf  LAPACK tridiagonal factorization routine.
35c     dgttrs  LAPACK tridiagonal solve routine.
36c     daxpy   Level 1 BLAS that computes y <- alpha*x+y.
37c     dcopy   Level 1 BLAS that copies one vector to another.
38c     dscal   Level 1 BLAS that scales a vector by a scalar.
39c     dnrm2   Level 1 BLAS that computes the norm of a vector. 
40c     av      Matrix vector multiplication routine that computes A*x.
41c     mv      Matrix vector multiplication routine that computes M*x.
42c
43c\Author
44c     Danny Sorensen
45c     Richard Lehoucq
46c     Chao Yang
47c     Dept. of Computational &
48c     Applied Mathematics
49c     Rice University 
50c     Houston, Texas 
51c 
52c\SCCS Information: @(#)
53c FILE: sdrv6.F   SID: 2.5   DATE OF SID: 10/17/00   RELEASE: 2
54c
55c\Remarks
56c     1. None
57c
58c\EndLib
59c------------------------------------------------------------------------
60c
61c     %-----------------------------%
62c     | Define leading dimensions   |
63c     | for all arrays.             |
64c     | MAXN:   Maximum dimension   |
65c     |         of the A allowed.   |
66c     | MAXNEV: Maximum NEV allowed |
67c     | MAXNCV: Maximum NCV allowed |
68c     %-----------------------------%
69c
70      integer          maxn, maxnev, maxncv, ldv
71      parameter        (maxn=256, maxnev=10, maxncv=25, 
72     &                 ldv=maxn)
73c
74c     %--------------%
75c     | Local Arrays |
76c     %--------------%
77c
78      Double precision
79     &                 v(ldv,maxncv), workl(maxncv*(maxncv+8)),
80     &                 workd(3*maxn), d(maxncv,2), resid(maxn), 
81     &                 ad(maxn), adl(maxn), adu(maxn), adu2(maxn),
82     &                 temp(maxn), ax(maxn), mx(maxn)
83      logical          select(maxncv)
84      integer          iparam(11), ipntr(11), ipiv(maxn)
85c
86c     %---------------%
87c     | Local Scalars |
88c     %---------------%
89c
90      character        bmat*1, which*2
91      integer          ido, n, nev, ncv, lworkl, info, j, ierr,
92     &                 nconv, maxitr, ishfts, mode
93      logical          rvec
94      Double precision
95     &                 sigma, r1, r2, tol, h
96c
97c     %------------%
98c     | Parameters |
99c     %------------%
100c
101      Double precision
102     &                 zero, one, two, four, six
103      parameter        (zero = 0.0D+0, one = 1.0D+0, 
104     &                  four = 4.0D+0, six = 6.0D+0,
105     &                  two = 2.0D+0 ) 
106c
107c     %-----------------------------%
108c     | BLAS & LAPACK routines used |
109c     %-----------------------------%
110c
111      Double precision
112     &                 dnrm2
113      external         daxpy, dcopy, dscal, dnrm2, dgttrf, dgttrs
114c
115c     %--------------------%
116c     | Intrinsic function |
117c     %--------------------%
118c
119      intrinsic        abs
120c
121c     %-----------------------%
122c     | Executable statements |
123c     %-----------------------%
124c
125c     %--------------------------------------------------%
126c     | The number N is the dimension of the matrix. A   |
127c     | generalized eigenvalue problem is solved (BMAT = |
128c     | 'G'.) NEV is the number of eigenvalues to be     |
129c     | approximated.  Since the Cayley mode is used,    |
130c     | WHICH is set to 'LM'.  The user can modify NEV,  |
131c     | NCV, SIGMA to solve problems of different sizes, |
132c     | and to get different parts of the spectrum.      |
133c     | However, The following conditions must be        |
134c     | satisfied:                                       |
135c     |                 N <= MAXN,                       | 
136c     |               NEV <= MAXNEV,                     |
137c     |           NEV + 1 <= NCV <= MAXNCV               | 
138c     %--------------------------------------------------%
139c
140      n = 100
141      nev = 4
142      ncv = 20
143      if ( n .gt. maxn ) then
144         print *, ' ERROR with _SDRV6: N is greater than MAXN '
145         go to 9000
146      else if ( nev .gt. maxnev ) then
147         print *, ' ERROR with _SDRV6: NEV is greater than MAXNEV '
148         go to 9000
149      else if ( ncv .gt. maxncv ) then
150         print *, ' ERROR with _SDRV6: NCV is greater than MAXNCV '
151         go to 9000
152      end if
153      bmat = 'G'
154      which = 'LM'
155      sigma = 1.5D+2
156c
157c     %--------------------------------------------------%
158c     | The work array WORKL is used in DSAUPD as        |
159c     | workspace.  Its dimension LWORKL is set as       |
160c     | illustrated below.  The parameter TOL determines |
161c     | the stopping criterion.  If TOL<=0, machine      |
162c     | precision is used.  The variable IDO is used for |
163c     | reverse communication and is initially set to 0. |
164c     | Setting INFO=0 indicates that a random vector is |
165c     | generated in DSAUPD to start the Arnoldi         |
166c     | iteration.                                       |
167c     %--------------------------------------------------%
168c
169      lworkl = ncv*(ncv+8)
170      tol = zero
171      ido = 0
172      info = 0
173c
174c     %---------------------------------------------------%
175c     | This program uses exact shifts with respect to    |
176c     | the current Hessenberg matrix (IPARAM(1) = 1).    |
177c     | IPARAM(3) specifies the maximum number of Arnoldi |
178c     | iterations allowed.  Mode 5 specified in the      |
179c     | documentation of DSAUPD is used (IPARAM(7) = 5).  |
180c     | All these options may be changed by the user. For |
181c     | details, see the documentation in DSAUPD.         |
182c     %---------------------------------------------------%
183c
184      ishfts = 1
185      maxitr = 300
186      mode   = 5
187c
188      iparam(1) = ishfts
189      iparam(3) = maxitr 
190      iparam(7) = mode 
191c
192c     %------------------------------------------------------%
193c     | Call LAPACK routine to factor (A-sigma*M).  The      |
194c     | stiffness matrix A is the 1-d discrete Laplacian.    |
195c     | The mass matrix M is the associated mass matrix      |
196c     | arising from using piecewise linear finite elements  |
197c     | on the interval [0, 1].                              |
198c     %------------------------------------------------------%
199c
200      h = one / dble(n+1)
201      r1 = (four / six) * h
202      r2 = (one / six) * h
203      do 20 j=1,n
204         ad(j) = two / h - sigma * r1
205         adl(j) = -one / h - sigma * r2
206 20   continue
207      call dcopy (n, adl, 1, adu, 1)
208      call dgttrf (n, adl, ad, adu, adu2, ipiv, ierr)
209      if (ierr .ne. 0) then
210         print *, ' '
211         print *, ' Error with _gttrf in _SDRV6.'
212         print *, ' '
213         go to 9000
214      end if
215c
216c     %-------------------------------------------%
217c     | M A I N   L O O P (Reverse communication) |
218c     %-------------------------------------------%
219c
220 10   continue
221c
222c        %---------------------------------------------%
223c        | Repeatedly call the routine DSAUPD and take |
224c        | actions indicated by parameter IDO until    |
225c        | either convergence is indicated or maxitr   |
226c        | has been exceeded.                          |
227c        %---------------------------------------------%
228c
229         call dsaupd ( ido, bmat, n, which, nev, tol, resid, ncv, 
230     &                 v, ldv, iparam, ipntr, workd, workl, lworkl, 
231     &                 info )
232c
233         if (ido .eq. -1) then
234c
235c           %-------------------------------------------------------%
236c           | Perform  y <--- OP*x = (inv[A-SIGMA*M])*(A+SIGMA*M)*x |
237c           | to force starting vector into the range of OP.  The   |
238c           | user should provide his/her matrix vector (A*x, M*x)  |
239c           | multiplication routines and a linear system solver    |
240c           | here.  The matrix vector multiplication routine takes |
241c           | workd(ipntr(1)) as the input vector.  The final       |
242c           | result is returned to workd(ipntr(2)).                | 
243c           %-------------------------------------------------------%
244c
245            call av (n, workd(ipntr(1)), workd(ipntr(2)))
246            call mv (n, workd(ipntr(1)), temp)
247            call daxpy(n, sigma, temp, 1, workd(ipntr(2)), 1) 
248c
249            call dgttrs ('Notranspose', n, 1, adl, ad, adu, adu2, ipiv, 
250     &                workd(ipntr(2)), n, ierr)
251            if (ierr .ne. 0) then
252               print *, ' '
253               print *, ' Error with _gttrs in _SDRV6.'
254               print *, ' ' 
255               go to 9000
256            end if
257c
258c           %-----------------------------------------%
259c           | L O O P   B A C K to call DSAUPD again. |
260c           %-----------------------------------------%
261c
262            go to 10
263c
264         else if (ido .eq. 1) then
265c
266c           %----------------------------------------------------%
267c           | Perform y <-- OP*x = inv[A-SIGMA*M]*(A+SIGMA*M)*x. |
268c           | M*x has been saved in workd(ipntr(3)).  The user   |
269c           | need the matrix vector multiplication (A*x)        |
270c           | routine and a linear system solver here.  The      |
271c           | matrix vector multiplication routine takes         |
272c           | workd(ipntr(1)) as the input, and the result is    |
273c           | combined with workd(ipntr(3)) to form the input    |
274c           | for the linear system solver.  The final result is |
275c           | returned to workd(ipntr(2)).                       | 
276c           %----------------------------------------------------%
277c
278            call av (n, workd(ipntr(1)), workd(ipntr(2)))
279            call daxpy(n, sigma, workd(ipntr(3)), 1, workd(ipntr(2)), 1)
280            call dgttrs ('Notranspose', n, 1, adl, ad, adu, adu2, ipiv, 
281     &                   workd(ipntr(2)), n, ierr)
282            if (ierr .ne. 0) then
283               print *, ' ' 
284               print *, ' Error with _gttrs in _SDRV6. '
285               print *, ' '
286               go to 9000
287            end if
288c
289c           %-----------------------------------------%
290c           | L O O P   B A C K to call DSAUPD again. |
291c           %-----------------------------------------%
292c
293            go to 10
294c
295         else if (ido .eq. 2) then
296c
297c           %--------------------------------------------%
298c           |             Perform  y <--- M*x.           |
299c           | Need matrix vector multiplication routine  |
300c           | here that takes workd(ipntr(1)) as input   |
301c           | and returns the result to workd(ipntr(2)). |
302c           %--------------------------------------------%
303c
304            call mv (n, workd(ipntr(1)), workd(ipntr(2)))
305c
306c           %-----------------------------------------%
307c           | L O O P   B A C K to call DSAUPD again. |
308c           %-----------------------------------------%
309c
310            go to 10
311c
312         end if
313c
314c     %-----------------------------------------%
315c     | Either we have convergence, or there is |
316c     | an error.                               |
317c     %-----------------------------------------%
318c
319      if ( info .lt. 0 ) then
320c
321c        %--------------------------%
322c        | Error message, check the |
323c        | documentation in DSAUPD  |
324c        %--------------------------%
325c
326           print *, ' '
327           print *, ' Error with _saupd, info = ',info
328           print *, ' Check the documentation of _saupd. '
329           print *, ' '
330c
331      else
332c
333c        %-------------------------------------------%
334c        | No fatal errors occurred.                 |
335c        | Post-Process using DSEUPD.                |
336c        |                                           |
337c        | Computed eigenvalues may be extracted.    |
338c        |                                           |
339c        | Eigenvectors may also be computed now if  |
340c        | desired.  (indicated by rvec = .true.)    |
341c        %-------------------------------------------%
342c
343         rvec = .true.
344c
345         call dseupd ( rvec, 'All', select, d, v, ldv, sigma, 
346     &        bmat, n, which, nev, tol, resid, ncv, v, ldv, 
347     &        iparam, ipntr, workd, workl, lworkl, ierr )
348c
349c        %----------------------------------------------%
350c        | Eigenvalues are returned in the first column |
351c        | of the two dimensional array D and the       |
352c        | corresponding eigenvectors are returned in   |
353c        | the first NEV columns of the two dimensional |
354c        | array V if requested.  Otherwise, an         |
355c        | orthogonal basis for the invariant subspace  |
356c        | corresponding to the eigenvalues in D is     |
357c        | returned in V.                               |
358c        %----------------------------------------------%
359c
360         if ( ierr .ne. 0 ) then
361c
362c           %------------------------------------%
363c           | Error condition:                   |
364c           | Check the documentation of DSEUPD. |
365c           %------------------------------------%
366c
367            print *, ' '
368            print *, ' Error with _seupd, info = ', ierr
369            print *, ' Check the documentation of _seupd '
370            print *, ' '
371c
372         else
373c
374c           %---------------------------%
375c           | Compute the residual norm |
376c           |                           |
377c           |   ||  A*x - lambda*x ||   |
378c           |                           |
379c           | for the NCONV accurately  |
380c           | computed eigenvalues and  |
381c           | eigenvectors.  (iparam(5) |
382c           | indicates how many are    |
383c           | accurate to the requested |
384c           | tolerance)                |
385c           %---------------------------%
386c 
387            nconv = iparam(5)
388            do 30 j=1, nconv
389               call av(n, v(1,j), ax)
390               call mv(n, v(1,j), mx)
391               call daxpy (n, -d(j,1), mx, 1, ax, 1)
392               d(j,2) = dnrm2(n, ax, 1)
393               d(j,2) = d(j,2) / abs(d(j,1))
394 30         continue
395c
396            call dmout(6, nconv, 2, d, maxncv, -6,
397     &           'Ritz values and relative residuals')
398c
399         end if
400c
401c        %------------------------------------------%
402c        | Print additional convergence information |
403c        %------------------------------------------%
404c
405         if ( info .eq. 1) then
406            print *, ' '
407            print *, ' Maximum number of iterations reached.'
408            print *, ' '
409         else if ( info .eq. 3) then
410            print *, ' '
411            print *, ' No shifts could be applied during implicit',
412     &               ' Arnoldi update, try increasing NCV.'
413            print *, ' '
414         end if
415c
416         print *, ' '
417         print *, ' _SDRV6 '
418         print *, ' ====== '
419         print *, ' '
420         print *, ' Size of the matrix is', n
421         print *, ' The number of Ritz values requested is', nev
422         print *, ' The number of Arnoldi vectors generated',
423     &            ' (NCV) is ', ncv
424         print *, ' What portion of the spectrum: ', which
425         print *, ' The number of converged Ritz values is ',
426     &              nconv 
427         print *, ' The number of Implicit Arnoldi update',
428     &           ' iterations taken is', iparam(3)
429         print *, ' The number of OP*x is ', iparam(9)
430         print *, ' The convergence criterion is ', tol
431         print *, ' '
432c
433      end if
434c
435c     %---------------------------%
436c     | Done with program dsdrv6. |
437c     %---------------------------%
438c
439 9000 continue
440c
441      end
442c
443c------------------------------------------------------------------------
444c     Matrix vector subroutine
445c     where the matrix used is the 1 dimensional mass matrix
446c     arising from using the piecewise linear finite element
447c     on the interval [0,1].
448c
449      subroutine mv (n, v, w)
450      integer           n, j
451      Double precision 
452     &                  v(n), w(n), one, four, six, h
453      parameter         (one = 1.0D+0, four = 4.0D+0, 
454     &                   six = 6.0D+0) 
455c
456      w(1) =  four*v(1) + v(2)
457      do 100 j = 2,n-1
458         w(j) = v(j-1) + four*v(j) + v(j+1) 
459  100 continue
460      j = n
461      w(j) = v(j-1) + four*v(j) 
462c
463c     Scale the vector w by h.
464c
465      h = one / (six*dble(n+1))
466      call dscal(n, h, w, 1)
467      return
468      end
469c
470c------------------------------------------------------------------------
471c     Matrix vector subroutine
472c     where the matrix is the stiffness matrix obtained from the 
473c     finite element discretization of the 1-dimensional discrete Laplacian
474c     on the interval [0,1] with zero Dirichlet boundary condition
475c     using piecewise linear elements.
476c
477      subroutine av (n, v, w)
478      integer           n, j
479      Double precision
480     &                  v(n), w(n), one, two, h
481      parameter         (one = 1.0D+0, two = 2.0D+0)
482c
483      w(1) =  two*v(1) - v(2)
484      do 100 j = 2,n-1
485         w(j) = - v(j-1) + two*v(j) - v(j+1) 
486  100 continue
487      j = n
488      w(j) = - v(j-1) + two*v(j)
489c     
490c     Scale the vector w by (1/h).
491c
492      h = one / dble(n+1)
493      call dscal(n, one/h, w, 1) 
494      return
495      end
Note: See TracBrowser for help on using the repository browser.