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