[2696] | 1 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2 | ! Driver to test Hessian Singular Vectors |
---|
| 3 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 4 | |
---|
| 5 | PROGRAM KPP_ROOT_Driver |
---|
| 6 | |
---|
| 7 | USE KPP_ROOT_Model |
---|
| 8 | USE KPP_ROOT_Initialize, ONLY: Initialize |
---|
| 9 | |
---|
| 10 | KPP_REAL :: T, DVAL(NSPEC), DV(NVAR) |
---|
| 11 | KPP_REAL :: Hs(NVAR,NVAR), SOA(NVAR) |
---|
| 12 | INTEGER :: i, j, ind_1 = 1, ind_2 = 2, ind_COST |
---|
| 13 | LOGICAL :: LHessian |
---|
| 14 | INTEGER, PARAMETER :: NSOA = 1 |
---|
| 15 | KPP_REAL :: Y_tlm(NVAR,NSOA) |
---|
| 16 | KPP_REAL :: Y_adj(NVAR,1) |
---|
| 17 | KPP_REAL :: Y_soa(NVAR,NSOA) |
---|
| 18 | |
---|
| 19 | integer lwork, n |
---|
| 20 | parameter (n = NVAR) |
---|
| 21 | |
---|
| 22 | ! Number of desired eigenvectors |
---|
| 23 | INTEGER, PARAMETER :: NEIG = 4 |
---|
| 24 | double complex zwork(2000), alpha(NEIG), beta(NEIG), & |
---|
| 25 | eivec(n,NEIG), tmp(n), residu(n) |
---|
| 26 | |
---|
| 27 | integer kmax, jmax, jmin, maxstep, method, m, l, maxnmv, & |
---|
| 28 | order, testspace |
---|
| 29 | double precision tol, lock, dznrm2 |
---|
| 30 | logical wanted |
---|
| 31 | double complex target |
---|
| 32 | real elapse |
---|
| 33 | real etime, tarray(2) |
---|
| 34 | |
---|
| 35 | target = (0.0,0.0) |
---|
| 36 | tol = 1.d-9 |
---|
| 37 | kmax = NEIG ! Number of wanted solutions |
---|
| 38 | jmin = 10; jmax = 20 ! min/max size of search space |
---|
| 39 | |
---|
| 40 | maxstep = 30 ! Max number of Jacobi-Davidson iterations |
---|
| 41 | lock = 1.d-9 |
---|
| 42 | !... order = 0: nearest to target |
---|
| 43 | !... order = -1: smallest real part |
---|
| 44 | !... order = 1: largest real part |
---|
| 45 | !... order = -2: smallest complex part |
---|
| 46 | !... order = 2: largest complex part |
---|
| 47 | order = 1 |
---|
| 48 | |
---|
| 49 | method = 2 ! 1=gmres(m), 2=cgstab(l) |
---|
| 50 | m = 30 ! for gmres(m): |
---|
| 51 | l = 2 ! for cgstab(l): |
---|
| 52 | |
---|
| 53 | maxnmv = 20 !... maximum number of matvecs in cgstab or gmres |
---|
| 54 | |
---|
| 55 | testspace = 3 |
---|
| 56 | !... Testspace 1: w = "Standard Petrov" * v (Section 3.1.1) |
---|
| 57 | !... Testspace 2: w = "Standard 'variable' Petrov" * v (Section 3.1.2) |
---|
| 58 | !... Testspace 3: w = "Harmonic Petrov" * v (Section 3.5.1) |
---|
| 59 | |
---|
| 60 | if ( method .eq. 1 ) then |
---|
| 61 | lwork = 4 + m + 5*jmax + 3*kmax |
---|
| 62 | else |
---|
| 63 | lwork = 10 + 6*l + 5*jmax + 3*kmax |
---|
| 64 | end if |
---|
| 65 | wanted = .true. |
---|
| 66 | |
---|
| 67 | call jdqz(alpha, beta, eivec, wanted, n, target, tol, & |
---|
| 68 | kmax, jmax, jmin, & |
---|
| 69 | method, m, l, maxnmv, maxstep, & |
---|
| 70 | lock, order, testspace, zwork, lwork ) |
---|
| 71 | |
---|
| 72 | elapse = etime( tarray ) |
---|
| 73 | |
---|
| 74 | print*,'Number of converged eigenvalues = ',kmax |
---|
| 75 | |
---|
| 76 | !... Compute the norms of the residuals: |
---|
| 77 | do j = 1, kmax |
---|
| 78 | call amul ( n, eivec(1,j), residu ) |
---|
| 79 | call zscal ( n, beta(j), residu, 1) |
---|
| 80 | call bmul ( n, eivec(1,j), tmp ) |
---|
| 81 | call zaxpy( n, -alpha(j), tmp, 1, residu, 1 ) |
---|
| 82 | print '("lambda(",i2,"): (",1p,e11.4,",",e11.4, & |
---|
| 83 | " )")', j,alpha(j)/beta(j) |
---|
| 84 | print '(a30,d13.6)', '||beta Ax - alpha Bx||:', & |
---|
| 85 | dznrm2( n, residu, 1 ) |
---|
| 86 | end do |
---|
| 87 | write(*,10) tarray(1), elapse |
---|
| 88 | |
---|
| 89 | 10 format(1x,'END JDQZ AFTER ',f6.2,' SEC. CPU-TIME AND ', f6.2, & |
---|
| 90 | ' SEC. ELAPSED TIME' ) |
---|
| 91 | |
---|
| 92 | open(120,file='KPP_ROOT_vec.m') |
---|
| 93 | do i=1,NVAR |
---|
| 94 | write(120,120) (eivec(i,j), j=1,NEIG) |
---|
| 95 | end do |
---|
| 96 | close(120) |
---|
| 97 | |
---|
| 98 | open(120,file='KPP_ROOT_val.m') |
---|
| 99 | do i=1,NEIG |
---|
| 100 | write(120,120) ALPHA(i), BETA(i), ALPHA(i)/BETA(i) |
---|
| 101 | end do |
---|
| 102 | close(120) |
---|
| 103 | |
---|
| 104 | 120 format(10000(E14.7,1X)) |
---|
| 105 | |
---|
| 106 | END PROGRAM KPP_ROOT_Driver |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | !========================================================================== |
---|
| 110 | ! Dummy preconditioner subroutine |
---|
| 111 | !========================================================================== |
---|
| 112 | subroutine PRECON( neq, q ) |
---|
| 113 | !............................................... |
---|
| 114 | !... Subroutine to compute q = K^-1 q |
---|
| 115 | !............................................... |
---|
| 116 | integer neq, i |
---|
| 117 | double complex q(neq) |
---|
| 118 | end subroutine PRECON |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | !========================================================================== |
---|
| 122 | ! matrix vector multiplication subroutine |
---|
| 123 | ! |
---|
| 124 | ! |
---|
| 125 | ! Compute the matrix vector multiplication y<---A*x |
---|
| 126 | ! where A = adjoint(Model)*tlm(Model) |
---|
| 127 | ! |
---|
| 128 | ! Note: we multiply the real and imaginaray parts separately |
---|
| 129 | ! x = (xr,xi) |
---|
| 130 | ! y = (yr,yi) |
---|
| 131 | ! We run the code twice such that: |
---|
| 132 | ! yr <-- A*xr and yi <-- A*xi |
---|
| 133 | ! |
---|
| 134 | !========================================================================== |
---|
| 135 | subroutine AMUL (n, u, v) |
---|
| 136 | USE KPP_ROOT_Model |
---|
| 137 | USE KPP_ROOT_Initialize, ONLY: Initialize |
---|
| 138 | integer n, j |
---|
| 139 | Double Complex u(n), v(n) |
---|
| 140 | Double precision xr(n), yr(n), xi(n), yi(n), three, two |
---|
| 141 | parameter (three = 3.0D+0, two = 2.0D+0) |
---|
| 142 | INTEGER, PARAMETER :: NSOA = 1 |
---|
| 143 | KPP_REAL :: Y_tlm(NVAR,NSOA) |
---|
| 144 | KPP_REAL :: Y_adj(NVAR,1) |
---|
| 145 | KPP_REAL :: Y_soa(NVAR,NSOA) |
---|
| 146 | integer, save :: indexA = 0 |
---|
| 147 | |
---|
| 148 | indexA = indexA + 1 |
---|
| 149 | print*,'AMUL #',indexA |
---|
| 150 | |
---|
| 151 | ! Real and imaginary parts of complex inputs |
---|
| 152 | xr(1:n) = DBLE( u(1:n) ) |
---|
| 153 | xi(1:n) = DIMAG( u(1:n) ) |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | DO i=1,NVAR |
---|
| 157 | RTOL(i) = 1.0d-4 |
---|
| 158 | ATOL(i) = 1.0d-3 |
---|
| 159 | END DO |
---|
| 160 | !~~~> The cost function is 0.5*VAR(ind_COST, tF)**2 |
---|
| 161 | ind_COST = ind_O3 |
---|
| 162 | |
---|
| 163 | ! First do a FWD and a TLM integration |
---|
| 164 | CALL Initialize() |
---|
| 165 | Y_tlm(1:NVAR,1) = xr(1:NVAR) |
---|
| 166 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 167 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 168 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 169 | ATOL, RTOL) |
---|
| 170 | |
---|
| 171 | ! Set Lambda(t_final) = Y_tlm(t_final) and do backward integration |
---|
| 172 | CALL Initialize() |
---|
| 173 | Y_adj(1:NVAR,1) = Y_tlm(1:NVAR) |
---|
| 174 | Y_tlm(1:NVAR,1) = 0.0d0 |
---|
| 175 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 176 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 177 | ATOL, RTOL) |
---|
| 178 | yr(1:NVAR) = Y_adj(1:NVAR) |
---|
| 179 | |
---|
| 180 | ! First do a FWD and a TLM integration |
---|
| 181 | CALL Initialize() |
---|
| 182 | Y_tlm(1:NVAR,1) = xi(1:NVAR) |
---|
| 183 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 184 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 185 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 186 | ATOL, RTOL) |
---|
| 187 | |
---|
| 188 | ! Set Lambda(t_final) = Y_tlm(t_final) and do backward integration |
---|
| 189 | CALL Initialize() |
---|
| 190 | Y_adj(1:NVAR,1) = Y_tlm(1:NVAR) |
---|
| 191 | Y_tlm(1:NVAR,1) = 0.0d0 |
---|
| 192 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 193 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 194 | ATOL, RTOL) |
---|
| 195 | yi(1:NVAR) = Y_adj(1:NVAR) |
---|
| 196 | |
---|
| 197 | ! Convert double outputs to complex |
---|
| 198 | DO i=1,n |
---|
| 199 | v(i) = DCMPLX( yr(i), yi(i) ) |
---|
| 200 | END DO |
---|
| 201 | |
---|
| 202 | |
---|
| 203 | print*,' (END AMUL #',indexA,')' |
---|
| 204 | |
---|
| 205 | end subroutine AMUL |
---|
| 206 | |
---|
| 207 | !========================================================================== |
---|
| 208 | ! matrix vector multiplication subroutine |
---|
| 209 | ! |
---|
| 210 | ! Compute the matrix vector multiplication y<---B*x |
---|
| 211 | ! where B = Hessian, and y = SOA |
---|
| 212 | ! |
---|
| 213 | ! Note: we multiply the real and imaginaray parts separately |
---|
| 214 | ! x = (xr,xi) |
---|
| 215 | ! y = (yr,yi) |
---|
| 216 | ! We run the code twice such that: |
---|
| 217 | ! yr <-- B*xr and yi <-- B*xi |
---|
| 218 | !========================================================================== |
---|
| 219 | |
---|
| 220 | subroutine BMUL (n, u, v) |
---|
| 221 | USE KPP_ROOT_Model |
---|
| 222 | USE KPP_ROOT_Initialize, ONLY: Initialize |
---|
| 223 | integer n, j |
---|
| 224 | Double Complex u(n), v(n) |
---|
| 225 | Double precision xr(n), yr(n), xi(n), yi(n) |
---|
| 226 | INTEGER, PARAMETER :: NSOA = 1 |
---|
| 227 | KPP_REAL :: Y_tlm(NVAR,NSOA) |
---|
| 228 | KPP_REAL :: Y_adj(NVAR,1) |
---|
| 229 | KPP_REAL :: Y_soa(NVAR,NSOA) |
---|
| 230 | KPP_REAL :: Final(NVAR), FinalTlm(NVAR) |
---|
| 231 | integer, save :: indexB = 0 |
---|
| 232 | |
---|
| 233 | indexB = indexB + 1 |
---|
| 234 | print*,'BMUL #',indexB |
---|
| 235 | |
---|
| 236 | ! Real and imaginary parts of complex inputs |
---|
| 237 | xr(1:n) = DBLE( u(1:n) ) |
---|
| 238 | xi(1:n) = DIMAG( u(1:n) ) |
---|
| 239 | |
---|
| 240 | DO i=1,NVAR |
---|
| 241 | RTOL(i) = 1.0d-4 |
---|
| 242 | ATOL(i) = 1.0d-3 |
---|
| 243 | END DO |
---|
| 244 | !~~~> The cost function is 0.5*VAR(ind_COST, tF)**2 |
---|
| 245 | ind_COST = ind_O3 |
---|
| 246 | |
---|
| 247 | ! Forward & TLM integration |
---|
| 248 | CALL Initialize() |
---|
| 249 | Y_tlm(1:NVAR,1) = xr |
---|
| 250 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 251 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 252 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 253 | ATOL, RTOL) |
---|
| 254 | Final(1:NVAR) = VAR(1:NVAR) |
---|
| 255 | FinalTlm(1:NVAR) = Y_tlm(1:NVAR) |
---|
| 256 | |
---|
| 257 | ! Adjoint and SOA integration, real part |
---|
| 258 | CALL Initialize() |
---|
| 259 | Y_adj(1:NVAR,1) = 0.0d0; Y_adj(ind_COST,1) = Final(ind_COST) |
---|
| 260 | Y_soa(1:NVAR,1) = 0.0d0; Y_soa(ind_COST,1) = FinalTlm(ind_COST) |
---|
| 261 | Y_tlm(1:NVAR,1) = xr |
---|
| 262 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 263 | ATOL, RTOL) |
---|
| 264 | |
---|
| 265 | yr(1:NVAR) = Y_soa(1:NVAR,1) |
---|
| 266 | |
---|
| 267 | |
---|
| 268 | ! Forward & TLM integration |
---|
| 269 | CALL Initialize() |
---|
| 270 | Y_tlm(1:NVAR,1) = xi |
---|
| 271 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 272 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 273 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 274 | ATOL, RTOL) |
---|
| 275 | Final(1:NVAR) = VAR(1:NVAR) |
---|
| 276 | FinalTlm(1:NVAR) = Y_tlm(1:NVAR) |
---|
| 277 | |
---|
| 278 | ! Adjoint and SOA integration, imaginary part |
---|
| 279 | CALL Initialize() |
---|
| 280 | Y_adj(1:NVAR,1) = 0.0d0; Y_adj(ind_COST,1) = Final(ind_COST) |
---|
| 281 | Y_soa(1:NVAR,1) = 0.0d0; Y_soa(ind_COST,1) = FinalTlm(ind_COST) |
---|
| 282 | Y_tlm(1:NVAR,1) = xi |
---|
| 283 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 284 | ATOL, RTOL) |
---|
| 285 | |
---|
| 286 | yi(1:NVAR) = Y_soa(1:NVAR,1) |
---|
| 287 | |
---|
| 288 | ! Convert double outputs to complex |
---|
| 289 | DO i=1,n |
---|
| 290 | v(i) = DCMPLX( yr(i), yi(i) ) |
---|
| 291 | END DO |
---|
| 292 | |
---|
| 293 | |
---|
| 294 | print*,' (END BMUL #',indexB,')' |
---|
| 295 | |
---|
| 296 | end subroutine BMUL |
---|
| 297 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 298 | |
---|
| 299 | |
---|