[2696] | 1 | program example |
---|
| 2 | c |
---|
| 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) |
---|
| 7 | c |
---|
| 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) |
---|
| 15 | c |
---|
| 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 |
---|
| 23 | c... order = 0: nearest to target |
---|
| 24 | c... order = -1: smallest real part |
---|
| 25 | c... order = 1: largest real part |
---|
| 26 | c... order = -2: smallest complex part |
---|
| 27 | c... order = 2: largest complex part |
---|
| 28 | order = 0 |
---|
| 29 | c... method = 1: gmres(m) |
---|
| 30 | c... method = 2: cgstab(l) |
---|
| 31 | method = 2 |
---|
| 32 | c... for gmres(m): |
---|
| 33 | m = 30 |
---|
| 34 | c... for cgstab(l): |
---|
| 35 | l= 2 |
---|
| 36 | c... maximum number of matvecs in cgstab or gmres |
---|
| 37 | maxnmv = 100 |
---|
| 38 | testspace = 3 |
---|
| 39 | c... Testspace 1: w = "Standard Petrov" * v (Section 3.1.1) |
---|
| 40 | c... Testspace 2: w = "Standard 'variable' Petrov" * v (Section 3.1.2) |
---|
| 41 | c... Testspace 3: w = "Harmonic Petrov" * v (Section 3.5.1) |
---|
| 42 | c |
---|
| 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. |
---|
| 49 | c |
---|
| 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 ) |
---|
| 54 | c |
---|
| 55 | elapse = etime( tarray ) |
---|
| 56 | c |
---|
| 57 | c... 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 |
---|
| 69 | c |
---|
| 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 ) |
---|
| 75 | c............................................... |
---|
| 76 | c... Subroutine to compute q = K^-1 q |
---|
| 77 | c............................................... |
---|
| 78 | integer neq, i |
---|
| 79 | double complex q(neq) |
---|
| 80 | c |
---|
| 81 | do i = 1, neq |
---|
| 82 | q(i) = i*q(i)/(i*i-31) |
---|
| 83 | end do |
---|
| 84 | c |
---|
| 85 | end |
---|
| 86 | c |
---|
| 87 | subroutine AMUL( neq, q, r ) |
---|
| 88 | c............................................... |
---|
| 89 | c... Subroutine to compute r = Aq |
---|
| 90 | c............................................... |
---|
| 91 | integer neq, i |
---|
| 92 | double complex q(neq), r(neq) |
---|
| 93 | c |
---|
| 94 | do i = 1, neq |
---|
| 95 | r(i) = i*q(i) |
---|
| 96 | end do |
---|
| 97 | c |
---|
| 98 | end |
---|
| 99 | c |
---|
| 100 | subroutine BMUL( neq, q, r ) |
---|
| 101 | c............................................... |
---|
| 102 | c... Subroutine to compute r = Bq |
---|
| 103 | c............................................... |
---|
| 104 | integer neq, i |
---|
| 105 | double complex q(neq), r(neq) |
---|
| 106 | c |
---|
| 107 | do i = 1, neq |
---|
| 108 | r(i) = q(i)/i |
---|
| 109 | end do |
---|
| 110 | c |
---|
| 111 | end |
---|
| 112 | |
---|