1 subroutine DEVVUN(A, LDA, N, EVALR, EVALI, VEC, IFLAG, WORK) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c This file contains DEVVUN and CDIV. 6c>> 1996-03-30 DEVUNN Krogh MIN0 => MIN, added external statement. 7c>> 1994-11-02 DEVVUN Krogh Changes to use M77CON 8c>> 1992-04-23 DEVVUN CLL Declaring all variables. 9c>> 1992-04-22 DEVVUN Krogh Removed unused labels 330 and 1220. 10c>> 1992-03-05 DEVVUN Krogh Initial version, converted from EISPACK. 11c 12c This subroutine uses slight modifcations of EISPACK routines 13c BALANC, ELMHES, ELTRAN, HQR2 and BALBAK to get the eigenvalues and 14c eigenvectors of a general real matrix by the QR method. The first 15c two are are encapsulated in DEVBH. The following three are given 16c inline here. 17c 18c ELTRAN is a translation of the algol procedure ELMTRANS, 19c Num. Math. 16, 181-204(1970) by Peters and Wilkinson. 20c Handbook for Auto. Comp., Vol.ii-Linear Algebra, 372-395(1971). 21c 22c HQR2 is a translation of the ALGOL procedure HQR2, 23c Num. Math. 16, 181-204(1970) by Peters and Wilkinson. 24c Handbook for Auto. Comp., Vol.ii-Linear Algebra, 372-395(1971). 25c 26c BALBAK is a translation of the ALGOL procedure BALBAK, 27c Num. Math. 13, 293-304(1969) by Parlett and Reinsch. 28c Handbook for Auto. Comp., Vol.ii-Linear Algebra, 315-326(1971). 29c 30c This subroutine finds the eigenvalues and eigenvectors 31c of a real general matrix by the QR method. 32c 33c On input 34c 35c LDA must be set to the row dimension of two-dimensional 36c array parameters as declared in the calling program 37c dimension statement. 38c N is the order of the matrix. 39c A contains the input matrix whose eigenvalues and eigenvectors 40c are desired. 41c 42c On output 43c A has been destroyed. 44c EVALR and EVALI contain real and imaginary parts, respectively, of 45c the eigenvalues. The eigenvalues are given in order of 46c increasing real parts. When real parts are equal they are 47c given in order of increasing absolute complex part. Complex 48c conjugate pairs of values appear consecutively with 49c the eigenvalue having the positive imaginary part first. If 50c an error exit is made, the eigenvalues should be correct 51c (but unordered) for indices IFLAG(1)+1,...,N. 52c VEC contains the real and imaginary parts of the eigenvectors. 53c if the I-th eigenvalue is real, the I-th column of VEC 54c contains its eigenvector. If the I-th eigenvalue is complex 55c with positive imaginary part, the I-th and (I+1)-th 56c columns of VEC contain the real and imaginary parts of its 57c eigenvector. The eigenvectors are unnormalized. If an 58c error exit is made, none of the eigenvectors has been found. 59c IFLAG(1) is set to 60c 1 If all eigenvalues are real. 61c 2 If some eigenvalues are complex. 62c 3 If N < 1 on the initial entry. 63c 4 If the limit of 30*N iterations is exhausted. 64c ------------------------------------------------------------------ 65c--D replaces "?": ?EVVUN, ?EVBH, ?NRM2, ?SCAL, ?CDIV 66c ------------------------------------------------------------------ 67 external DNRM2 68 integer I,J,K,L,M,N,EN,NA,LDA,IGH,ITN,ITS,LOW,MP,MP2,ENM2 69 integer II, IFLAG(N), LTYPE 70 double precision A(LDA,N),EVALR(N),EVALI(N),VEC(LDA,N), WORK(N) 71 double precision P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2 72 double precision DNRM2 73 logical NOTLAS 74c ------------------------------------------------------------------ 75c 76 call DEVBH(A, LDA, N, LOW, IGH, IFLAG, WORK) 77 LTYPE = 1 78c 79c -------------------------------- ELTRAN --------------------------- 80c 81c Accumulate the stabilized elementary similarity transformations 82c used in the preceding reduction of the matrix to upper Hessenberg 83c form. 84c 85c .......... initialize VEC to identity matrix .......... 86 do 20 J = 1, N 87 do 10 I = 1, N 88 VEC(I,J) = 0.0D0 89 10 continue 90 VEC(J,J) = 1.0D0 91 20 continue 92 if (IGH .gt. LOW+1) then 93 do 50 MP = IGH-1, LOW+1, -1 94 do 30 I = MP+1, IGH 95 VEC(I,MP) = A(I,MP-1) 96 30 continue 97 I = IFLAG(MP) 98 if (I .ne. MP) then 99 do 40 J = MP, IGH 100 VEC(MP,J) = VEC(I,J) 101 VEC(I,J) = 0.0D0 102 40 continue 103 VEC(I,MP) = 1.0D0 104 end if 105 50 continue 106 end if 107c 108c -------------------------------- HQR2 ----------------------------- 109c 110c Find the eigenvalues of a real upper Hessenberg matrix by the QR 111c method. 112c 113 IFLAG(1) = 0 114 NORM = 0.0D0 115 K = 1 116c .......... store roots isolated by balanc 117c and compute matrix norm .......... 118 do 70 I = 1, N 119 do 60 J = K, N 120 NORM = NORM + abs(A(I,J)) 121 60 continue 122 K = I 123 if ((I .lt. LOW) .or. (I .gt. IGH)) then 124 EVALR(I) = A(I,I) 125 EVALI(I) = 0.0D0 126 end if 127 70 continue 128 EN = IGH 129 T = 0.0D0 130 ITN = 30*N 131c .......... search for next eigenvalues .......... 132 80 if (EN .lt. LOW) go to 340 133 ITS = 0 134 NA = EN - 1 135 ENM2 = NA - 1 136c .......... look for single small sub-diagonal element 137 90 do 100 L = EN, LOW+1, -1 138 S = abs(A(L-1,L-1)) + abs(A(L,L)) 139 if (S .eq. 0.0D0) S = NORM 140 TST1 = S 141 TST2 = TST1 + abs(A(L,L-1)) 142 if (TST2 .eq. TST1) go to 110 143 100 continue 144c .......... form shift .......... 145 110 X = A(EN,EN) 146 if (L .eq. EN) go to 270 147 Y = A(NA,NA) 148 W = A(EN,NA) * A(NA,EN) 149 if (L .eq. NA) go to 280 150 if (ITN .le. 0) then 151c .......... set error -- all eigenvalues have not 152c converged after 30*N iterations .......... 153 call ERMSG('DEVVUN', EN, 0, 154 1'ERROR NO. is index of eigenvalue causing convergence failure.', 155 2 '.') 156 IFLAG(1) = 4 157 if (EN .le. 0) IFLAG(1) = 3 158 return 159 end if 160 if (ITS .ne. 10 .and. ITS .ne. 20) go to 130 161c .......... form exceptional shift .......... 162 T = T + X 163 do 120 I = LOW, EN 164 A(I,I) = A(I,I) - X 165 120 continue 166 S = abs(A(EN,NA)) + abs(A(NA,ENM2)) 167 X = 0.75D0 * S 168 Y = X 169 W = -0.4375D0 * S * S 170 130 ITS = ITS + 1 171 ITN = ITN - 1 172c .......... look for two consecutive small sub-diagonal elements. 173 do 140 M = ENM2, L, -1 174 ZZ = A(M,M) 175 R = X - ZZ 176 S = Y - ZZ 177 P = (R * S - W) / A(M+1,M) + A(M,M+1) 178 Q = A(M+1,M+1) - ZZ - R - S 179 R = A(M+2,M+1) 180 S = abs(P) + abs(Q) + abs(R) 181 P = P / S 182 Q = Q / S 183 R = R / S 184 if (M .eq. L) go to 150 185 TST1 = abs(P)*(abs(A(M-1,M-1)) + abs(ZZ) + abs(A(M+1,M+1))) 186 TST2 = TST1 + abs(A(M,M-1))*(abs(Q) + abs(R)) 187 if (TST2 .eq. TST1) go to 150 188 140 continue 189 150 MP2 = M + 2 190 do 160 I = MP2, EN 191 A(I,I-2) = 0.0D0 192 if (I .ne. MP2) A(I,I-3) = 0.0D0 193 160 continue 194c .......... double QR step involving rows L to EN and 195c columns M to EN .......... 196 do 260 K = M, NA 197 NOTLAS = K .ne. NA 198 if (K .ne. M) then 199 P = A(K,K-1) 200 Q = A(K+1,K-1) 201 R = 0.0D0 202 if (NOTLAS) R = A(K+2,K-1) 203 X = abs(P) + abs(Q) + abs(R) 204 if (X .eq. 0.0D0) go to 260 205 P = P / X 206 Q = Q / X 207 R = R / X 208 end if 209 S = sign(sqrt(P*P+Q*Q+R*R),P) 210 if (K .ne. M) then 211 A(K,K-1) = -S * X 212 else 213 if (L .ne. M) A(K,K-1) = -A(K,K-1) 214 end if 215 P = P + S 216 X = P / S 217 Y = Q / S 218 ZZ = R / S 219 Q = Q / P 220 R = R / P 221 if (NOTLAS) then 222c .......... row modification .......... 223 do 200 J = K, N 224 P = A(K,J) + Q * A(K+1,J) + R * A(K+2,J) 225 A(K,J) = A(K,J) - P * X 226 A(K+1,J) = A(K+1,J) - P * Y 227 A(K+2,J) = A(K+2,J) - P * ZZ 228 200 continue 229 J = min(EN,K+3) 230c .......... column modification .......... 231 do 210 I = 1, J 232 P = X * A(I,K) + Y * A(I,K+1) + ZZ * A(I,K+2) 233 A(I,K) = A(I,K) - P 234 A(I,K+1) = A(I,K+1) - P * Q 235 A(I,K+2) = A(I,K+2) - P * R 236 210 continue 237c .......... accumulate transformations .......... 238 do 220 I = LOW, IGH 239 P = X * VEC(I,K) + Y * VEC(I,K+1) + ZZ * VEC(I,K+2) 240 VEC(I,K) = VEC(I,K) - P 241 VEC(I,K+1) = VEC(I,K+1) - P * Q 242 VEC(I,K+2) = VEC(I,K+2) - P * R 243 220 continue 244 else 245c .......... row modification .......... 246 do 230 J = K, N 247 P = A(K,J) + Q * A(K+1,J) 248 A(K,J) = A(K,J) - P * X 249 A(K+1,J) = A(K+1,J) - P * Y 250 230 continue 251c 252 J = min(EN,K+3) 253c .......... column modification .......... 254 do 240 I = 1, J 255 P = X * A(I,K) + Y * A(I,K+1) 256 A(I,K) = A(I,K) - P 257 A(I,K+1) = A(I,K+1) - P * Q 258 240 continue 259c .......... accumulate transformations .......... 260 do 250 I = LOW, IGH 261 P = X * VEC(I,K) + Y * VEC(I,K+1) 262 VEC(I,K) = VEC(I,K) - P 263 VEC(I,K+1) = VEC(I,K+1) - P * Q 264 250 continue 265 end if 266 260 continue 267 go to 90 268c .......... one root found .......... 269 270 A(EN,EN) = X + T 270 EVALR(EN) = A(EN,EN) 271 EVALI(EN) = 0.0D0 272 EN = NA 273 go to 80 274c .......... two roots found .......... 275 280 P = (Y - X) / 2.0D0 276 Q = P * P + W 277 ZZ = sqrt(abs(Q)) 278 A(EN,EN) = X + T 279 X = A(EN,EN) 280 A(NA,NA) = Y + T 281 if (Q .ge. 0.0D0) then 282c .......... real pair .......... 283 ZZ = P + sign(ZZ,P) 284 EVALR(NA) = X + ZZ 285 EVALR(EN) = EVALR(NA) 286 if (ZZ .ne. 0.0D0) EVALR(EN) = X - W / ZZ 287 EVALI(NA) = 0.0D0 288 EVALI(EN) = 0.0D0 289 X = A(EN,NA) 290 S = abs(X) + abs(ZZ) 291 P = X / S 292 Q = ZZ / S 293 R = sqrt(P*P+Q*Q) 294 P = P / R 295 Q = Q / R 296c .......... row modification .......... 297 do 290 J = NA, N 298 ZZ = A(NA,J) 299 A(NA,J) = Q * ZZ + P * A(EN,J) 300 A(EN,J) = Q * A(EN,J) - P * ZZ 301 290 continue 302c .......... column modification .......... 303 do 300 I = 1, EN 304 ZZ = A(I,NA) 305 A(I,NA) = Q * ZZ + P * A(I,EN) 306 A(I,EN) = Q * A(I,EN) - P * ZZ 307 300 continue 308c .......... accumulate transformations .......... 309 do 310 I = LOW, IGH 310 ZZ = VEC(I,NA) 311 VEC(I,NA) = Q * ZZ + P * VEC(I,EN) 312 VEC(I,EN) = Q * VEC(I,EN) - P * ZZ 313 310 continue 314 else 315c .......... complex pair .......... 316 LTYPE = 2 317 EVALR(NA) = X + P 318 EVALR(EN) = X + P 319 EVALI(NA) = ZZ 320 EVALI(EN) = -ZZ 321 end if 322 EN = ENM2 323 go to 80 324c .......... all roots found. Backsubstitute to find 325c vectors of upper triangular form .......... 326 340 if (NORM .eq. 0.0D0) go to 1000 327 do 800 EN = N, 1, -1 328 P = EVALR(EN) 329 Q = EVALI(EN) 330 NA = EN - 1 331 if (Q) 710, 600, 800 332c .......... real vector .......... 333 600 M = EN 334 A(EN,EN) = 1.0D0 335 do 700 I = NA, 1, -1 336 W = A(I,I) - P 337 R = 0.0D0 338 do 610 J = M, EN 339 R = R + A(I,J) * A(J,EN) 340 610 continue 341 if (EVALI(I) .lt. 0.0D0) then 342 ZZ = W 343 S = R 344 else 345 M = I 346 if (EVALI(I) .eq. 0.0D0) then 347 T = W 348 if (T .eq. 0.0D0) then 349 TST1 = NORM 350 T = TST1 351 640 T = 0.01D0 * T 352 TST2 = NORM + T 353 if (TST2 .gt. TST1) go to 640 354 end if 355 A(I,EN) = -R / T 356 else 357c .......... solve real equations .......... 358 X = A(I,I+1) 359 Y = A(I+1,I) 360 Q = (EVALR(I)-P) * (EVALR(I)-P) + EVALI(I) * EVALI(I) 361 T = (X * S - ZZ * R) / Q 362 A(I,EN) = T 363 if (abs(X) .gt. abs(ZZ)) then 364 A(I+1,EN) = (-R - W * T) / X 365 else 366 A(I+1,EN) = (-S - Y * T) / ZZ 367 end if 368 end if 369c .......... overflow control .......... 370 T = abs(A(I,EN)) 371 if (T .ne. 0.0D0) then 372 TST1 = T 373 TST2 = TST1 + 1.0D0/TST1 374 if (TST2 .le. TST1) then 375 do 690 J = I, EN 376 A(J,EN) = A(J,EN)/T 377 690 continue 378 end if 379 end if 380 end if 381 700 continue 382c .......... end real vector .......... 383 go to 800 384c 385c .......... complex vector .......... 386 710 M = NA 387c .......... last vector component chosen imaginary so that 388c eigenvector matrix is triangular .......... 389 if (abs(A(EN,NA)) .gt. abs(A(NA,EN))) then 390 A(NA,NA) = Q / A(EN,NA) 391 A(NA,EN) = -(A(EN,EN) - P) / A(EN,NA) 392 else 393 call DCDIV(0.0D0,-A(NA,EN),A(NA,NA)-P,Q,A(NA,NA),A(NA,EN)) 394 end if 395 A(EN,NA) = 0.0D0 396 A(EN,EN) = 1.0D0 397 ENM2 = NA - 1 398 do 760 I = ENM2, 1, -1 399 W = A(I,I) - P 400 RA = 0.0D0 401 SA = 0.0D0 402 do 730 J = M, EN 403 RA = RA + A(I,J) * A(J,NA) 404 SA = SA + A(I,J) * A(J,EN) 405 730 continue 406 if (EVALI(I) .lt. 0.0D0) then 407 ZZ = W 408 R = RA 409 S = SA 410 else 411 M = I 412 if (EVALI(I) .eq. 0.0D0) then 413 call DCDIV(-RA,-SA,W,Q,A(I,NA),A(I,EN)) 414 else 415c .......... solve complex equations .......... 416 X = A(I,I+1) 417 Y = A(I+1,I) 418 VR = (EVALR(I)-P)*(EVALR(I)-P)+EVALI(I)*EVALI(I)-Q*Q 419 VI = (EVALR(I) - P) * 2.0D0 * Q 420 if (VR .eq. 0.0D0 .and. VI .eq. 0.0D0) then 421 TST1 = NORM * (abs(W)+abs(Q)+abs(X)+abs(Y)+abs(ZZ)) 422 VR = TST1 423 740 VR = 0.01D0 * VR 424 TST2 = TST1 + VR 425 if (TST2 .gt. TST1) go to 740 426 end if 427 call DCDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI, 428 x A(I,NA),A(I,EN)) 429 if (abs(X) .gt. abs(ZZ) + abs(Q)) then 430 A(I+1,NA) = (-RA - W * A(I,NA) + Q * A(I,EN)) / X 431 A(I+1,EN) = (-SA - W * A(I,EN) - Q * A(I,NA)) / X 432 else 433 call DCDIV(-R-Y*A(I,NA),-S-Y*A(I,EN),ZZ,Q, 434 x A(I+1,NA),A(I+1,EN)) 435 end if 436 end if 437c .......... overflow control .......... 438 T = max(abs(A(I,NA)), abs(A(I,EN))) 439 if (T .ne. 0.0D0) then 440 TST1 = T 441 TST2 = TST1 + 1.0D0/TST1 442 if (TST2 .le. TST1) then 443 do 750 J = I, EN 444 A(J,NA) = A(J,NA)/T 445 A(J,EN) = A(J,EN)/T 446 750 continue 447 end if 448 end if 449 end if 450 760 continue 451c .......... end complex vector .......... 452 800 continue 453c .......... end back substitution. 454c vectors of isolated roots .......... 455 do 840 I = 1, N 456 if (I .ge. LOW .and. I .le. IGH) go to 840 457c 458 do 820 J = I, N 459 VEC(I,J) = A(I,J) 460 820 continue 461c 462 840 continue 463c .......... multiply by transformation matrix to give 464c vectors of original full matrix. 465 do 880 J = N, LOW, -1 466 M = min(J,IGH) 467c 468 do 880 I = LOW, IGH 469 ZZ = 0.0D0 470c 471 do 860 K = LOW, M 472 ZZ = ZZ + VEC(I,K) * A(K,J) 473 860 continue 474c 475 VEC(I,J) = ZZ 476 880 continue 477c 478 1000 continue 479c 480c ------------------------------ BALBAK ----------------------------- 481c 482c Form eigenvectors of a real general matrix by back transforming 483c those of the corresponding balanced matrix determined by BALANC. 484c 485 if (IGH .ne. LOW) then 486 do 1110 I = LOW, IGH 487 S = WORK(I) 488c .......... left hand eigenvectors are back transformed 489c if the foregoing statement is replaced by 490c S=1.0D0/WORK(I). .......... 491 do 1100 J = 1, N 492 VEC(I,J) = VEC(I,J) * S 493 1100 continue 494 1110 continue 495 end if 496c ......... for I=LOW-1 step -1 until 1, 497c IGH+1 step 1 until N do -- .......... 498 do 1140 II = 1, N 499 I = II 500 if ((I .lt. LOW) .or. (I .gt. IGH)) then 501 if (I .lt. LOW) I = LOW - II 502 K = WORK(I) 503 if (K .ne. I) then 504 do 1130 J = 1, N 505 S = VEC(I,J) 506 VEC(I,J) = VEC(K,J) 507 VEC(K,J) = S 508 1130 continue 509 end if 510 end if 511 1140 continue 512c Normalize the eigenvectors 513 do 1220 I = 1, N 514 P = DNRM2(N, VEC(1, I), 1) 515 if (EVALI(I) .eq. 0.D0) then 516 call DSCAL(N, sign(1.D0, VEC(1,I)) / P, VEC(1,I), 1) 517 else if (EVALI(I) .gt. 0.D0) then 518 Q = DNRM2(N, VEC(1, I+1), 1) 519 if (P .gt. Q) then 520 P = P * sqrt(1.D0 + (Q/P)**2) 521 else if (Q .ne. 0.D0) then 522 P = Q * sqrt(1.D0 + (P/Q)**2) 523 else 524 go to 1220 525 end if 526 if (VEC(1, I+1) .eq. 0.D0) then 527 P = sign(1.D0, VEC(1,I+1)) / P 528 Q = 0.D0 529 else if (abs(VEC(1, I)) .gt. abs(VEC(1, I+1))) then 530 P = 1.D0 / (P*sqrt(1.D0+(VEC(1,I+1)/VEC(1,I))**2)) 531 Q = P * (VEC(1,I+1) / abs(VEC(1,I))) 532 P = sign(P, VEC(1, I)) 533 else 534 Q = 1.D0 / (P*sqrt(1.D0+(VEC(1,I)/VEC(1,I+1))**2)) 535 P = Q * (VEC(1,I) / abs(VEC(1,I+1))) 536 Q = sign(Q, VEC(1, I+1)) 537 end if 538 do 1200 J = 1, N 539 R = P * VEC(J, I) + Q * VEC(J, I+1) 540 VEC(J, I+1) = P * VEC(J, I+1) - Q * VEC(J, I) 541 VEC(J, I) = R 542 1200 continue 543 end if 544 1220 continue 545c-- Begin mask code changes 546c Set up for Shell sort 547c Sort so real parts are algebraically increasing 548c For = real parts, so abs(imag. parts) are increasing 549c For both =, sort on index -- preserves complex pair order 550c-- End mask code changes 551 do 2000 I = 1, N 552 IFLAG(I) = I 553 2000 continue 554 L = 1 555 do 2010 K = 1, N 556 L = 3*L + 1 557 if (L .ge. N) go to 2020 558 2010 continue 559 2020 L = max(1, (L-4) / 9) 560 2030 do 2100 J = L+1, N 561 K = IFLAG(J) 562 P = EVALR(K) 563 I = J - L 564 2040 if (P - EVALR(IFLAG(I))) 2070, 2050, 2080 565 2050 if (abs(EVALI(K)) - abs(EVALI(IFLAG(I)))) 2070, 2060, 2080 566 2060 if (K .gt. IFLAG(I)) go to 2080 567 2070 IFLAG(I+L) = IFLAG(I) 568 I = I - L 569 if (I .gt. 0) go to 2040 570 2080 IFLAG(I+L) = K 571 2100 continue 572 L = (L-1) / 3 573 if (L .ne. 0) go to 2030 574c Indices in IFLAG now give the desired order -- 575c Move entries to get this order. 576 2110 do 2150 I = L+1, N 577 if (IFLAG(I) .ne. I) then 578 L = I 579 M = I 580 P = EVALR(I) 581 Q = EVALI(I) 582 do 2115 J = 1, N 583 WORK(J) = VEC(J, I) 584 2115 continue 585 2120 K = IFLAG(M) 586 IFLAG(M) = M 587 if (K .ne. L) then 588 EVALR(M) = EVALR(K) 589 EVALI(M) = EVALI(K) 590 do 2130 J = 1, N 591 VEC(J, M) = VEC(J, K) 592 2130 continue 593 M = K 594 go to 2120 595 else 596 EVALR(M) = P 597 EVALI(M) = Q 598 do 2140 J = 1, N 599 VEC(J, M) = WORK(J) 600 2140 continue 601 go to 2110 602 end if 603 end if 604 2150 continue 605 IFLAG(1) = LTYPE 606 return 607 end 608c ================================================================== 609 subroutine DCDIV(AR,AI,BR,BI,CR,CI) 610c 611c complex division, (CR,CI) = (AR,AI)/(BR,BI) 612c ------------------------------------------------------------------ 613 double precision AR,AI,BR,BI,CR,CI 614 double precision S,ARS,AIS,BRS,BIS 615 S = abs(BR) + abs(BI) 616 ARS = AR/S 617 AIS = AI/S 618 BRS = BR/S 619 BIS = BI/S 620 S = BRS**2 + BIS**2 621 CR = (ARS*BRS + AIS*BIS)/S 622 CI = (AIS*BRS - ARS*BIS)/S 623 return 624 end 625