source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/examples/example.f @ 4239

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

Merge of branch palm4u into trunk

File size: 3.1 KB
Line 
1      program example
2c
3      integer          lwork, n
4      parameter        (n = 100)
5      double complex   zwork(200 000), alpha(20), beta(20), 
6     $                 eivec(n,20), tmp(n), residu(n)
7c
8      integer          kmax, jmax, jmin, maxstep, method, m, l, maxnmv, 
9     $                 order, testspace, j
10      double precision tol, lock, dznrm2
11      logical          wanted
12      double complex   target
13      real             elapse
14      real             etime, tarray(2)
15c
16      target = (31.,0.0)
17      tol = 1.d-9
18      kmax = 5
19      jmax = 20
20      jmin = 10
21      maxstep = 1000
22      lock = 1.d-9
23c...     order =  0: nearest to target
24c...     order = -1: smallest real part
25c...     order =  1: largest real part
26c...     order = -2: smallest complex part
27c...     order =  2: largest complex part
28      order = 0
29c...     method = 1: gmres(m)
30c...     method = 2: cgstab(l)
31      method = 2
32c...     for gmres(m):
33      m = 30
34c...     for cgstab(l):
35      l= 2
36c...     maximum number of matvecs in cgstab or gmres
37      maxnmv = 100
38      testspace = 3
39c...     Testspace 1: w = "Standard Petrov" * v (Section 3.1.1)
40c...     Testspace 2: w = "Standard 'variable' Petrov" * v (Section 3.1.2)
41c...     Testspace 3: w = "Harmonic Petrov" * v (Section 3.5.1)
42c
43      if ( method .eq. 1 ) then
44         lwork =  4 +  m  + 5*jmax + 3*kmax
45      else
46         lwork = 10 + 6*l + 5*jmax + 3*kmax
47      end if
48      wanted = .true.
49c
50      call jdqz(alpha, beta, eivec, wanted, n, target, tol, 
51     $     kmax, jmax, jmin,
52     $     method, m, l, maxnmv, maxstep,
53     $     lock, order, testspace, zwork, lwork )
54c
55      elapse = etime( tarray )
56c
57c...     Compute the norms of the residuals:
58      do j = 1, kmax
59         call amul  ( n, eivec(1,j), residu )
60         call zscal ( n, beta(j), residu, 1)
61         call bmul  ( n, eivec(1,j), tmp )
62         call zaxpy( n, -alpha(j), tmp, 1, residu, 1 )
63         print '("lambda(",i2,"): (",1p,e11.4,",",e11.4,
64     $        " )")', j,alpha(j)/beta(j)
65         print '(a30,d13.6)', '||beta Ax - alpha Bx||:',
66     $          dznrm2( n, residu, 1 )
67      end do
68      write(*,10) tarray(1), elapse
69c
70   10 format(1x,'END JDQZ AFTER ',f6.2,' SEC. CPU-TIME AND ', f6.2,
71     $       ' SEC. ELAPSED TIME' )
72      end
73
74      subroutine PRECON( neq, q )
75c...............................................
76c...     Subroutine to compute q = K^-1 q
77c...............................................
78      integer        neq, i
79      double complex q(neq)
80c
81      do i = 1, neq
82         q(i) = i*q(i)/(i*i-31)
83      end do
84c
85      end
86c
87      subroutine AMUL( neq, q, r )
88c...............................................
89c...     Subroutine to compute r = Aq
90c...............................................
91      integer        neq, i
92      double complex q(neq), r(neq)
93c
94      do i = 1, neq
95         r(i) = i*q(i)
96      end do
97c
98      end
99c
100      subroutine BMUL( neq, q, r )
101c...............................................
102c...     Subroutine to compute r = Bq
103c...............................................
104      integer        neq, i
105      double complex q(neq), r(neq)
106c
107      do i = 1, neq
108         r(i) = q(i)/i
109      end do
110c
111      end
112
Note: See TracBrowser for help on using the repository browser.