1C 2C modlib.f: grab-bag of mathematical functions 3C 4C This library is free software: you can redistribute it and/or 5C modify it under the terms of the GNU Lesser General Public License 6C version 3, modified in accordance with the provisions of the 7C license to address the requirements of UK law. 8C 9C You should have received a copy of the modified GNU Lesser General 10C Public License along with this library. If not, copies may be 11C downloaded from http://www.ccp4.ac.uk/ccp4license.php 12C 13C This program is distributed in the hope that it will be useful, 14C but WITHOUT ANY WARRANTY; without even the implied warranty of 15C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16C GNU Lesser General Public License for more details. 17C 18C cross.f dot.f ea06c.f ea08c.f ea09c.f 19C fa01as.f fa01bs.f fa01cs.f fa01ds.f fm02ad.f 20C icross.f idot.f iminv3 match.f matmul.f 21C matmulnm.f matmulgen.f 22C matvec.f mc04b.f minvn.f minv3.f ranmar.f 23C scalev.f transp.f unit.f vdif.f vset.f 24C vsum.f zipin.f zipout.f 25C 26C Routines added by pjx (June 2000): 27C 28C GMPRD multiply two matrices (general) 29C MATMUL4 multiply two 4x4 matrices 30C MATMLN multiply two nxn matrices 31C DETMAT calculate the determinant of 3x3 matrix 32C MATVC4 matrix times vector (matrix is 4x4 array, treated as 3x3, 33C vector is length 3) 34C ML3MAT multiplies three matrices of any size 35C INV44 invert 4x4 matrix 36C MATMLI integer version of 3x3 matrix product MATMUL 37C MATMULTRANS multiply 3x3 matrix by transpose of another 3x3 matrix 38C IMATVEC integer version of MATVEC (post-multiply 3x3 matrix with a 39C vector) 40C TRANSFRM variant of MATVC4, same except that the input vector is 41C overwritten by the output vector 42C 43C Routines added by mdw (April 2003) 44C 45C ZJVX Compute zero of Bessel function of 1st kind Jv(x) 46C JDVX Compute Bessel functions of 1st kind Jv(x) and their derivatives 47C JVX Compute Bessel functions of 1st kind Jv(x) 48C GAMMA Compute gamma function GA(x) 49C MSTA1 Determine the starting point for backward 50C recurrence such that the magnitude of 51C Jn(x) at that point is about 10^(-MP) 52C MSTA2 Determine the starting point for backward 53C recurrence such that all Jn(x) has MP 54C significant digits 55C 56C_BEGIN_GMPRD 57C 58 SUBROUTINE GMPRD(A,B,R,N,M,L) 59C ============================= 60C 61C 62C---- Ssp general matrix product 63C 64C R(N,L) = A(N,M) * B(M,L) 65C 66C .. Scalar Arguments .. 67 INTEGER L,M,N 68C .. 69C .. Array Arguments .. 70 REAL A(N*M),B(M*L),R(N*L) 71C .. 72C .. Local Scalars .. 73 INTEGER I,IB,IK,IR,J,JI,K 74C .. 75C 76C_END_GMPRD 77C 78 IR = 0 79 IK = -M 80 DO 30 K = 1,L 81 IK = IK + M 82 DO 20 J = 1,N 83 IR = IR + 1 84 JI = J - N 85 IB = IK 86 R(IR) = 0.0 87 DO 10 I = 1,M 88 JI = JI + N 89 IB = IB + 1 90 R(IR) = A(JI)*B(IB) + R(IR) 91 10 CONTINUE 92 20 CONTINUE 93 30 CONTINUE 94 END 95C 96C_BEGIN_MATMUL4 97C 98C ========================= 99 SUBROUTINE MATMUL4(A,B,C) 100C ========================= 101C 102C Multiply two 4x4 matrices, A=B*C 103C 104C 105 IMPLICIT NONE 106C 107C .. Array Arguments .. 108 REAL A(4,4),B(4,4),C(4,4) 109C .. 110C .. Local Scalars .. 111 REAL S 112 INTEGER I,J,K 113C .. 114C_END_MATMUL4 115C 116 DO 30 I = 1,4 117 DO 20 J = 1,4 118 S = 0 119 DO 10 K = 1,4 120 S = B(I,K)*C(K,J) + S 121 10 CONTINUE 122 A(I,J) = S 123 20 CONTINUE 124 30 CONTINUE 125 RETURN 126 END 127C 128C_BEGIN_MATMLN 129C 130 subroutine matmln(n,a,b,c) 131C ========================= 132C 133C Multiply two nxn matrices 134C a = b . c 135C 136 integer n 137 real a(n,n),b(n,n),c(n,n) 138 integer i,j,k 139C 140C_END_MATMLN 141C 142 do 1 i=1,n 143 do 2 j=1,n 144 a(j,i)=0. 145 do 3 k=1,n 146 a(j,i)= b(j,k)*c(k,i)+a(j,i) 147 3 continue 148 2 continue 149 1 continue 150 return 151 end 152C 153C_BEGIN_DETMAT 154C 155C =========================== 156 SUBROUTINE DETMAT(RMAT,DET) 157C ============================ 158C 159C---- Calculate determinant DET of 3x3 matrix RMAT 160C 161C 162C .. Scalar Arguments .. 163 REAL DET 164C .. 165C .. Array Arguments .. 166 REAL RMAT(3,3) 167C .. 168C_END_DETMAT 169C 170 DET = RMAT(1,1)*RMAT(2,2)*RMAT(3,3) + 171 + RMAT(1,2)*RMAT(2,3)*RMAT(3,1) + 172 + RMAT(1,3)*RMAT(2,1)*RMAT(3,2) - 173 + RMAT(1,3)*RMAT(2,2)*RMAT(3,1) - 174 + RMAT(1,2)*RMAT(2,1)*RMAT(3,3) - 175 + RMAT(1,1)*RMAT(2,3)*RMAT(3,2) 176C 177 END 178C 179C 180C_BEGIN_MATVC4 181C 182 SUBROUTINE MATVC4(A,R,B) 183C ======================== 184C 185C Matrix x Vector 186C A(3) = R(4x4) . B(3) 187C 188 REAL A(3), R(4,4), B(3) 189C 190 REAL AA(4), BB(4) 191C 192 EXTERNAL GMPRD,VSET 193C 194C_END_MATVC4 195C 196 CALL VSET(BB,B) 197 BB(4) = 1.0 198 CALL GMPRD(R,BB,AA,4,4,1) 199 CALL VSET(A,AA) 200 RETURN 201 END 202C 203C 204C_BEGIN_ML3MAT 205C 206C 26-Nov-1988 J. W. Pflugrath Cold Spring Harbor Laboratory 207C Edited to conform to Fortran 77. Renamed from Multiply_3_matrices to 208C ML3MAT 209C 210C ============================================================================== 211C 212C! to multiply three matrices 213C 214 SUBROUTINE ML3MAT 215C ! input: 1st side of 1st matrix 216 1 (P 217C ! input: first matrix 218 2 ,A 219C ! input: 2nd side of 1st matrix & 1st side of 2nd matrix 220 3 ,Q 221C ! input: second matrix 222 4 ,B 223C ! input: 2nd side of 2nd matrix & 1st side of 3rd matrix 224 5 ,R 225C ! input: third matrix 226 6 ,C 227C ! input: 2nd side of 3rd matrix 228 7 ,S 229C ! output: product matrix 230 8 ,D) 231C 232CEE Multiplies three real matrices of any dimensions. It is not optimised 233C for very large matrices. 234C Multiply_3_matrices 235C*** this routine is inefficient! 236C Multiply_3_matrices 237Created: 15-NOV-1985 D.J.Thomas, MRC Laboratory of Molecular Biology, 238C Hills Road, Cambridge, CB2 2QH, England 239C 240C ! loop counters 241 INTEGER I,J,K,L 242C ! loop limits 243 INTEGER P,Q,R,S 244C ! first input matrix 245 REAL A (1:P,1:Q) 246C ! second input matrix 247 REAL B (1:Q,1:R) 248C ! third input matrix 249 REAL C (1:R,1:S) 250C ! output matrix 251 REAL D (1:P,1:S) 252C 253C_END_ML3MAT 254C 255 DO 100 L = 1, S 256 DO 100 I = 1, P 257 D(I,L) = 0.0 258 DO 100 K = 1, R 259 DO 100 J = 1, Q 260C 261C ! accumulate product matrix D=ABC 262C 263100 D(I,L) = D(I,L) + A(I,J) * B(J,K) * C(K,L) 264 CONTINUE 265 CONTINUE 266 CONTINUE 267 CONTINUE 268C Multiply_3_matrices 269 RETURN 270 END 271C 272C 273C_BEGIN_INV44 274C 275C ====================== 276 SUBROUTINE INV44(A,AI) 277C ====================== 278C 279C SUBROUTINE TO INVERT 4*4 MATRICES FOR CONVERSION BETWEEN 280C FRACTIONAL AND ORTHOGONAL AXES 281C PARAMETERS 282C 283C A (I) 4*4 MATRIX TO BE INVERTED 284C AI (O) INVERSE MATRIX 285C 286C SPECIFICATION STATEMENTS 287C ------------------------ 288C 289C 290C GET COFACTORS OF 'A' IN ARRAY 'C' 291C --------------------------------- 292C 293 IMPLICIT NONE 294C 295C .. Array Arguments .. 296 REAL A(4,4),AI(4,4) 297C .. 298C .. Local Scalars .. 299 REAL AM,D 300 INTEGER I,I1,II,J,J1,JJ 301C .. 302C .. Local Arrays .. 303 REAL C(4,4),X(3,3) 304C .. 305C 306C_END_INV44 307C 308 DO 40 II = 1,4 309 DO 30 JJ = 1,4 310 I = 0 311 DO 20 I1 = 1,4 312 IF (I1.EQ.II) GO TO 20 313 I = I + 1 314 J = 0 315 DO 10 J1 = 1,4 316 IF (J1.EQ.JJ) GO TO 10 317 J = J + 1 318 X(I,J) = A(I1,J1) 319 10 CONTINUE 320 20 CONTINUE 321 AM = X(1,1)*X(2,2)*X(3,3) - X(1,1)*X(2,3)*X(3,2) + 322 + X(1,2)*X(2,3)*X(3,1) - X(1,2)*X(2,1)*X(3,3) + 323 + X(1,3)*X(2,1)*X(3,2) - X(1,3)*X(2,2)*X(3,1) 324 C(II,JJ) = (-1)** (II+JJ)*AM 325 30 CONTINUE 326 40 CONTINUE 327C 328C---- calculate determinant 329C 330 D = 0 331 DO 50 I = 1,4 332 D = A(I,1)*C(I,1) + D 333 50 CONTINUE 334C 335C---- get inverse matrix 336C 337 DO 70 I = 1,4 338 DO 60 J = 1,4 339 AI(I,J) = C(J,I)/D 340 60 CONTINUE 341 70 CONTINUE 342 RETURN 343 END 344C 345C 346C_BEGIN_MATMLI 347C 348 SUBROUTINE MATMLI(A,B,C) 349C ======================== 350C Integer matrix multiply 351 INTEGER A(3,3),B(3,3),C(3,3),I,J,K 352C 353C_END_MATMULI 354C 355 DO 1 I=1,3 356 DO 2 J=1,3 357 A(I,J)=0 358 DO 3 K=1,3 359 A(I,J)=A(I,J)+B(I,K)*C(K,J) 3603 CONTINUE 361 2 CONTINUE 362 1 CONTINUE 363 RETURN 364 END 365C 366C 367C_BEGIN_MATMULTRANS 368C 369 SUBROUTINE MATMULTrans(A,B,C) 370C ============================= 371C 372C A=B*C(transpose) for 3x3 matrices 373C 374 IMPLICIT NONE 375C 376C ..Array arguments 377 REAL A(3,3),B(3,3),C(3,3) 378C ..Local arrays 379 REAL CT(3,3) 380C 381 EXTERNAL MATMUL 382C 383C_END_MATMULTRANS 384C 385 CALL TRANSP(CT,C) 386 CALL MATMUL(A,B,CT) 387 RETURN 388 END 389C 390C 391C_BEGIN_IMATVEC 392C 393 SUBROUTINE IMATVEC(V,A,B) 394C ======================== 395C 396C---- Post-multiply a 3x3 matrix by a vector 397C 398C V=AB 399C 400C .. Array Arguments .. 401 INTEGER A(3,3),B(3),V(3) 402C .. 403C .. Local Scalars .. 404 INTEGER I,J,S 405C 406C_END_IMATVEC 407C 408 DO 20 I = 1,3 409 S = 0 410C 411C 412 DO 10 J = 1,3 413 S = A(I,J)*B(J) + S 414 10 CONTINUE 415C 416C 417 V(I) = S 418 20 CONTINUE 419C 420C 421 END 422C 423C 424C_BEGIN_TRANSFRM 425C 426 SUBROUTINE TRANSFRM(X,MAT) 427C ========================== 428C 429C Transform vector X(3) by quine matrix MAT(4,4) 430C Return transformed vector in X. 431C 432 IMPLICIT NONE 433C 434C ..Array arguments.. 435 REAL X(3),MAT(4,4) 436C ..Local arrays.. 437 REAL TMP(3) 438C 439C_END_TRANSFRM 440C 441 CALL MATVC4(TMP,MAT,X) 442 CALL VSET(X,TMP) 443 RETURN 444 END 445C 446C 447C 448C_BEGIN_CROSS 449C 450 SUBROUTINE CROSS(A,B,C) 451C ======================= 452C 453C compute vector product A = B x C 454C 455C .. Array Arguments .. 456 REAL A(3),B(3),C(3) 457C 458C_END_CROSS 459C .. 460 A(1) = B(2)*C(3) - C(2)*B(3) 461 A(2) = B(3)*C(1) - C(3)*B(1) 462 A(3) = B(1)*C(2) - C(1)*B(2) 463 END 464C 465C 466C_BEGIN_DOT 467C 468 REAL FUNCTION DOT(A,B) 469C ====================== 470C 471C dot product of two vectors 472C 473C .. Array Arguments .. 474 REAL A(3),B(3) 475C 476C_END_DOT 477C .. 478 DOT = A(1)*B(1) + A(2)*B(2) + A(3)*B(3) 479 END 480C 481 482C ****************************************************************** 483 SUBROUTINE EIGEN_RS_ASC(A, R, N, MV) 484C ****************************************************************** 485C 486C---- SUBROUTINE TO COMPUTE EIGENVALUES & EIGENVECTORS OF A REAL 487C---- SYMMETRIC MATRIX, FROM IBM SSP MANUAL (SEE P165). 488C---- DESCRIPTION OF PARAMETERS - 489C---- A - ORIGINAL MATRIX STORED COLUMNWISE AS UPPER TRIANGLE ONLY, 490C---- I.E. "STORAGE MODE" = 1. EIGENVALUES ARE WRITTEN INTO DIAGONAL 491C---- ELEMENTS OF A I.E. A(1) A(3) A(6) FOR A 3*3 MATRIX. 492C---- R - RESULTANT MATRIX OF EIGENVECTORS STORED COLUMNWISE IN SAME 493C---- ORDER AS EIGENVALUES. 494C---- N - ORDER OF MATRICES A & R. 495C---- MV = 0 TO COMPUTE EIGENVALUES & EIGENVECTORS. 496C 497C 498c IMPLICIT NONE 499 500 INTEGER NMAX 501 PARAMETER (NMAX=10) 502 503 REAL A(*), R(*) 504 INTEGER N, MV 505C 506 INTEGER IQ, J, I, IJ, IA, IND, L, M, MQ, LQ, LM, LL, MM 507 INTEGER ILQ, IMQ, IM, IL, ILR, IMR 508 REAL ANORM, ANRMX, RANGE, THR, X, Y, SINX, SINX2, 509 & COSX, COSX2, SINCS 510 511C Lapack variables 512 LOGICAL LAPACK 513 CHARACTER*1 JOBZ, LAPRANGE, UPLO 514 INTEGER INFO, NVECTORS 515 REAL ABSTOL 516 INTEGER ISUPPZ(2*NMAX), IWORK(10*NMAX) 517 REAL WORK(26*NMAX),EVALUES(NMAX),AM(NMAX,NMAX) 518C 519C-- FOR REAL 520C DATA RANGE/1D-12/ 521 DATA RANGE/1E-6/ 522 523C Alternative lapack routine - only marginally tested 524 LAPACK = .FALSE. 525 IF (LAPACK) THEN 526 NVECTORS = 0 527 IF (N.GT.NMAX) 528 + CALL CCPERR(1,'s/r EIGEN_RS_ASC: redimension NMAX!') 529 IF (MV.EQ.0) JOBZ = 'V' 530 LAPRANGE = 'A' 531 UPLO = 'U' 532 ABSTOL = 0.0 533 IA = 0 534 DO I = 1,N 535 DO J = 1,I 536 IA = IA + 1 537 AM(J,I) = A(IA) 538 ENDDO 539 ENDDO 540C CALL SSYEVR(JOBZ, LAPRANGE, UPLO, N, AM, N, 541C + 1, N, 1, N, ABSTOL, NVECTORS, EVALUES, R, 542C + N, ISUPPZ, WORK, 26*N, IWORK, 10*N, INFO) 543 IA = 0 544 DO I = 1,NVECTORS 545 DO J = 1,I 546 IA = IA + 1 547 IF (J.EQ.I) A(IA) = EVALUES(I) 548 ENDDO 549 ENDDO 550 RETURN 551 ENDIF 552 553 IF (MV.EQ.0) THEN 554 IQ=-N 555 DO J=1,N 556 IQ=IQ+N 557 DO I=1,N 558 IJ=IQ+I 559 IF (I.EQ.J) THEN 560 R(IJ)=1. 561 ELSE 562 R(IJ)=0. 563 ENDIF 564 ENDDO 565 ENDDO 566 ENDIF 567C 568C---- INITIAL AND FINAL NORMS (ANORM & ANRMX) 569 IA=0 570 ANORM=0. 571 DO I=1,N 572 DO J=1,I 573 IA=IA+1 574 IF (J.NE.I) ANORM=ANORM+A(IA)**2 575 ENDDO 576 ENDDO 577C 578 IF (ANORM.LE.0.) GOTO 165 579 ANORM=SQRT(2.*ANORM) 580 ANRMX=ANORM*RANGE/N 581C 582C---- INITIALIZE INDICATORS AND COMPUTE THRESHOLD 583 IND=0 584 THR=ANORM 58545 THR=THR/N 58650 L=1 58755 LQ=L*(L-1)/2 588 LL=L+LQ 589 M=L+1 590 ILQ=N*(L-1) 591C 592C---- COMPUTE SIN & COS 59360 MQ=M*(M-1)/2 594 LM=L+MQ 595 IF (A(LM)*A(LM)-THR.LT.0.) GOTO 130 596 IND=1 597 MM=M+MQ 598 X=.5*(A(LL)-A(MM)) 599 Y=-A(LM)/SQRT(A(LM)**2+X*X) 600C---- Protect against rounding error. C. Flensburg 20080307. 601 IF (ABS(Y).GT.1.0) Y=SIGN(1.0,Y) 602 IF (X.LT.0.) Y=-Y 603 SINX=Y/SQRT(2.*(1.+(SQRT(1.-Y*Y)))) 604 SINX2=SINX**2 605 COSX=SQRT(1.-SINX2) 606 COSX2=COSX**2 607 SINCS=SINX*COSX 608C 609C---- ROTATE L & M COLUMNS 610 IMQ=N*(M-1) 611 DO 125 I=1,N 612 IQ=I*(I-1)/2 613 IF (I.NE.L .AND. I.NE.M) THEN 614 IF (I.LT.M) THEN 615 IM=I+MQ 616 ELSE 617 IM=M+IQ 618 ENDIF 619 IF (I.LT.L) THEN 620 IL=I+LQ 621 ELSE 622 IL=L+IQ 623 ENDIF 624 X=A(IL)*COSX-A(IM)*SINX 625 A(IM)=A(IL)*SINX+A(IM)*COSX 626 A(IL)=X 627 ENDIF 628 IF (MV.EQ.0) THEN 629 ILR=ILQ+I 630 IMR=IMQ+I 631 X=R(ILR)*COSX-R(IMR)*SINX 632 R(IMR)=R(ILR)*SINX+R(IMR)*COSX 633 R(ILR)=X 634 ENDIF 635125 CONTINUE 636C 637 X=2.*A(LM)*SINCS 638 Y=A(LL)*COSX2+A(MM)*SINX2-X 639 X=A(LL)*SINX2+A(MM)*COSX2+X 640 A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) 641 A(LL)=Y 642 A(MM)=X 643C 644C---- TESTS FOR COMPLETION 645C---- TEST FOR M = LAST COLUMN 646130 IF (M.NE.N) THEN 647 M=M+1 648 GOTO 60 649 ENDIF 650C 651C---- TEST FOR L =PENULTIMATE COLUMN 652 IF (L.NE.N-1) THEN 653 L=L+1 654 GOTO55 655 ENDIF 656 IF (IND.EQ.1) THEN 657 IND=0 658 GOTO50 659 ENDIF 660C 661C---- COMPARE THRESHOLD WITH FINAL NORM 662 IF (THR.GT.ANRMX) GOTO 45 663165 RETURN 664 END 665C 666C The routines ea06c, ea08c, ea09c, fa01as, fa01bs, fa01cs, fa01ds, 667C fm02ad, mc04b, (and possibly others) are from the 668C Harwell Subroutine library. The conditions on their external use, 669C reproduced from the Harwell manual are: 670C * due acknowledgement is made of the use of subroutines in any 671C research publications resulting from their use. 672C * the subroutines may be modified for use in research applications 673C by external users. The nature of such modifiactions should be 674C indicated in writing for information to the liaison officer. At 675C no time however, shall the subroutines or modifications thereof 676C become the property of the external user. 677C The liaison officer for the library's external affairs is listed 678C as: Mr. S. Marlow, Building 8.9, Harwell Laboratory, Didcot, 679C Oxon OX11 0RA, UK. 680C 681C_BEGIN_EA06C 682C 683 SUBROUTINE EA06C(A,VALUE,VECTOR,M,IA,IV,W) 684C ========================================== 685C 686C** 18/03/70 LAST LIBRARY UPDATE 687C 688C ( Calls EA08C(W,W(M1),VALUE,VECTOR,M,IV,W(M+M1)) 689C and MC04B(A,W,W(M1),M,IA,W(M+M1)) ) 690C 691C Given a real MxM symmetric matrix A = {aij} this routine 692C finds all its eigenvalues (lambda)i i=1,2,.....,m and 693C eigenvectors xj i=1,2,...,m. i.e. finds the non-trivial 694C solutions of Ax=(lambda)x 695C 696C The matrix is reduced to tri-diagonal form by applying 697C Householder transformations. The eigenvalue problem for 698C the reduced problem is then solved using the QR algorithm 699C by calling EA08C. 700C 701C Argument list 702C ------------- 703C 704C IA (I) (integer) should be set to the first dimension 705C of the array A, i.e. if the allocation 706C for the array A was specified by 707C DIMENSION A(100,50) 708C then IA would be set to 100 709C 710C M (I) (integer) should be set to the order m of the matrix 711C 712C IV (I) (integer) should be set to the first dimension 713C of the 2-dimensional array VECTOR 714C 715C VECTOR(IV,M) (O) (real) 2-dimensional array, with first dimension IV, 716C containing the eigenvectors. The components 717C of the eigenvector vector(i) corresponding 718C to the eigenvalue (lambda)i (in VALUE(I)) 719C are placed in VECTOR(J,I) J=1,2,...,M. 720C The eigenvectors are normalized so that 721C xT(i)x(i)=1 i=1,2,...,m. 722C 723C VALUE(M) (O) (real) array in which the routine puts 724C the eigenvalues (lambda)i, i=1,2,...,m. 725C These are not necessarily in any order. 726C 727C W (I) (real(*)) working array used by the routine for 728C work space. The dimension must be set 729C to at least 5*M. 730C 731C_END_EA06C 732C 733C 734C .. Scalar Arguments .. 735 INTEGER IA,IV,M 736C .. 737C .. Array Arguments .. 738 REAL A(IA,M),VALUE(M),VECTOR(IV,M),W(*) 739C .. 740C .. Local Scalars .. 741 REAL PP 742 INTEGER I,I1,II,K,L,M1 743C .. 744C .. External Subroutines .. 745 EXTERNAL EA08C,MC04B 746C .. 747 M1 = M + 1 748 W(1) = A(1,1) 749 IF (M-2) 30,10,20 750 10 W(2) = A(2,2) 751 W(4) = A(2,1) 752 GO TO 30 753 754 20 CALL MC04B(A,W,W(M1),M,IA,W(M+M1)) 755 30 CALL EA08C(W,W(M1),VALUE,VECTOR,M,IV,W(M+M1)) 756 IF (M.GT.2) THEN 757 DO 70 L = 1,M 758 DO 60 II = 3,M 759 I = M - II + 1 760 IF (W(M1+I).NE.0) THEN 761 PP = 0.0 762 I1 = I + 1 763 DO 40 K = I1,M 764 PP = A(I,K)*VECTOR(K,L) + PP 765 40 CONTINUE 766 PP = PP/ (A(I,I+1)*W(M1+I)) 767 DO 50 K = I1,M 768 VECTOR(K,L) = A(I,K)*PP + VECTOR(K,L) 769 50 CONTINUE 770 END IF 771 772 60 CONTINUE 773 70 CONTINUE 774 END IF 775 776 END 777C 778C_BEGIN_EA08C 779C 780 SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W) 781C ===================================== 782C 783C (Calls EA09C(A,B,W(M+1),M,W)) 784C 785C This uses QR iteration to find the all the eigenvalues and 786C eigenvectors of the real symmetric tri-diagonal matrix 787C whose diagonal elements are A(i), i=1,M and off-diagonal 788C elements are B(i),i=2,M. The eigenvalues will have unit 789C length. The array W is used for workspace and must have 790C dimension at least 2*M. We treat VEC as if it had 791C dimensions (IV,M). 792C 793C First EA09, which uses the QR algorithm, is used to find 794C the eigenvalues; using these as shifts the QR algorithm is 795C again applied but now using the plane rotations to generate 796C the eigenvectors. Finally the eigenvalues are refined 797C by taking Rayleigh quotients of the vectors. 798C 799C Argument list 800C ------------- 801C 802C A(M) (I) (real) Diagonal elements 803C 804C B(M) (I) (real) Off-diagonal elements 805C 806C IV (I) (integer) should be set to the first dimen- 807C sion of the 2-dimensional array VEC 808C 809C M (I) (integer) should be set to the order m of the 810C matrix 811C 812C VALUE(M) (O) (real) Eigenvalues 813C 814C VEC (O) (real) Eigenvectors. The dimensions 815C should be set to (IV,M). 816C 817C W(*) (I) (real) Working array.The dimension must be 818C set to at least 2*M. 819C 820C_END_EA08C 821C 822C .. Scalar Arguments .. 823 INTEGER IV,M 824C .. 825C .. Array Arguments .. 826 REAL A(M),B(M),VALUE(M),VEC(1),W(*) 827C .. 828C .. Local Scalars .. 829 REAL A11,A12,A13,A21,A22,A23,A33,A34,CO,EPS,ROOT,S,SI, 830 + V1,V2,XAX,XX 831 INTEGER I,II,ITER,J,J1,J2,JI,JK,K,L,N1,N2,N2M1 832C .. 833C .. External Subroutines .. 834 EXTERNAL EA09C 835C .. 836C .. Intrinsic Functions .. 837 INTRINSIC ABS,MIN,SIGN,SQRT 838C .. 839C .. Data statements .. 840 DATA EPS/1.0E-5/,A34/0.0/ 841C .. 842C 843C 844 CALL EA09C(A,B,W(M+1),M,W) 845C 846C---- Set vec to the identity matrix. 847C 848 DO 20 I = 1,M 849 VALUE(I) = A(I) 850 W(I) = B(I) 851 K = (I-1)*IV + 1 852 L = K + M - 1 853 DO 10 J = K,L 854 VEC(J) = 0.0 855 10 CONTINUE 856 VEC(K+I-1) = 1.0 857 20 CONTINUE 858 ITER = 0 859 IF (M.NE.1) THEN 860 N2 = M 861 30 CONTINUE 862C 863C---- Each qr iteration is performed of rows and columns n1 to n2 864C 865 DO 40 II = 2,N2 866 N1 = 2 + N2 - II 867 IF (ABS(W(N1)).LE. (ABS(VALUE(N1-1))+ABS(VALUE(N1)))* 868 + EPS) GO TO 50 869 40 CONTINUE 870 N1 = 1 871 50 IF (N2.NE.N1) THEN 872 ROOT = W(M+N2) 873 ITER = ITER + 1 874 N2M1 = N2 - 1 875 A22 = VALUE(N1) 876 A12 = A22 - ROOT 877 A23 = W(N1+1) 878 A13 = A23 879 DO 70 I = N1,N2M1 880 A33 = VALUE(I+1) 881 IF (I.NE.N2M1) A34 = W(I+2) 882 S = SIGN(SQRT(A12*A12+A13*A13),A12) 883 SI = A13/S 884 CO = A12/S 885 JK = I*IV + 1 886 J1 = JK - IV 887 J2 = MIN(M,I+ITER) + J1 - 1 888 DO 60 JI = J1,J2 889 V1 = VEC(JI) 890 V2 = VEC(JK) 891 VEC(JI) = V1*CO + V2*SI 892 VEC(JK) = V2*CO - V1*SI 893 JK = JK + 1 894 60 CONTINUE 895 IF (I.NE.N1) W(I) = S 896 A11 = CO*A22 + SI*A23 897 A12 = CO*A23 + SI*A33 898 A13 = SI*A34 899 A21 = CO*A23 - SI*A22 900 A22 = CO*A33 - SI*A23 901 A23 = CO*A34 902 VALUE(I) = A11*CO + A12*SI 903 A12 = -A11*SI + A12*CO 904 W(I+1) = A12 905 A22 = A22*CO - A21*SI 906 70 CONTINUE 907 VALUE(N2) = A22 908 GO TO 30 909 910 ELSE 911 N2 = N2 - 1 912 IF (N2-1.GT.0) GO TO 30 913 END IF 914C 915C---- Rayleigh quotient 916C 917 DO 90 J = 1,M 918 K = (J-1)*IV 919 XX = VEC(K+1)**2 920 XAX = A(1)*XX 921 DO 80 I = 2,M 922 XX = VEC(K+I)**2 + XX 923 XAX = (B(I)*2.0*VEC(K+I-1)+VEC(K+I)*A(I))*VEC(K+I) + 924 + XAX 925 80 CONTINUE 926 VALUE(J) = XAX/XX 927 90 CONTINUE 928 END IF 929 930 END 931C 932C 933 SUBROUTINE EA09C(A,B,VALUE,M,OFF) 934C ================================= 935C 936C 18/03/70 LAST LIBRARY UPDATE 937C 938C .. Scalar Arguments .. 939 INTEGER M 940C .. 941C .. Array Arguments .. 942 REAL A(M),B(M),OFF(M),VALUE(M) 943C .. 944C .. Local Scalars .. 945 REAL A11,A12,A13,A21,A22,A23,A33,A34,BB,CC,CO,EPS, 946 + ROOT,S,SBB,SI 947 INTEGER I,II,N1,N2,N2M1 948C .. 949C .. Intrinsic Functions .. 950 INTRINSIC ABS,SQRT 951C .. 952C .. Data statements .. 953 DATA A34/0.0/,EPS/0.6E-7/ 954C .. 955C 956C 957 VALUE(1) = A(1) 958 IF (M.NE.1) THEN 959 DO 10 I = 2,M 960 VALUE(I) = A(I) 961 OFF(I) = B(I) 962 10 CONTINUE 963C 964C---- Each qr iteration is performed of rows and columns n1 to n2 965C 966 N2 = M 967 20 CONTINUE 968 IF (N2.GT.1) THEN 969 DO 30 II = 2,N2 970 N1 = 2 + N2 - II 971 IF (ABS(OFF(N1)).LE. (ABS(VALUE(N1-1))+ 972 + ABS(VALUE(N1)))*EPS) GO TO 40 973 30 CONTINUE 974 N1 = 1 975 40 IF (N2.NE.N1) THEN 976C 977C---- Root is the eigenvalue of the bottom 2*2 matrix that is nearest 978C to the last matrix element and is used to accelerate the 979C convergence 980C 981 BB = (VALUE(N2)-VALUE(N2-1))*0.5 982 CC = OFF(N2)*OFF(N2) 983 SBB = 1.0 984 IF (BB.LT.0.0) SBB = -1.0 985 ROOT = CC/ (SQRT(BB*BB+CC)*SBB+BB) + VALUE(N2) 986 N2M1 = N2 - 1 987 A22 = VALUE(N1) 988 A12 = A22 - ROOT 989 A23 = OFF(N1+1) 990 A13 = A23 991 DO 50 I = N1,N2M1 992 A33 = VALUE(I+1) 993 IF (I.NE.N2M1) A34 = OFF(I+2) 994 S = SQRT(A12*A12+A13*A13) 995 SI = A13/S 996 CO = A12/S 997 IF (I.NE.N1) OFF(I) = S 998 A11 = CO*A22 + SI*A23 999 A12 = CO*A23 + SI*A33 1000 A13 = SI*A34 1001 A21 = CO*A23 - SI*A22 1002 A22 = CO*A33 - SI*A23 1003 A23 = CO*A34 1004 VALUE(I) = A11*CO + A12*SI 1005 A12 = -A11*SI + A12*CO 1006 OFF(I+1) = A12 1007 A22 = A22*CO - A21*SI 1008 50 CONTINUE 1009 VALUE(N2) = A22 1010 1011 ELSE 1012 N2 = N2 - 1 1013 END IF 1014 1015 GO TO 20 1016 1017 END IF 1018 1019 END IF 1020 1021 END 1022C 1023C 1024C---- Changes put in to make it work on vax (this version has scaling) 1025C 1. 12 statements for calculation of over/underflow of determinant in ma21 1026C replaced by a simple one which will overflow more easily(entry ma21cd) 1027C 2. changes to mc10ad to replace 370-specific parts.... 1028C a. u=floati (6 times) 1029C b. new alog16 procedure (twice) 1030C c. simpler 16**diag statement (once) 1031C 3. all double precision replaced by real*8 (not really necessary). 1032C 4. replace a(n),b(n) by a(1),b(1) in fm02ad to avoid vax array checking. 1033C 1034 FUNCTION FA01AS(I) 1035C ================== 1036 DOUBLE PRECISION G, FA01AS 1037 INTEGER I 1038 COMMON/FA01ES/G 1039 G= 1431655765.D0 1040C 1041C 1042 G=DMOD(G* 9228907.D0,4294967296.D0) 1043 IF(I.GE.0)FA01AS=G/4294967296.D0 1044 IF(I.LT.0)FA01AS=2.D0*G/4294967296.D0-1.D0 1045 RETURN 1046 END 1047C 1048 SUBROUTINE FA01BS(MAX,NRAND) 1049C ============================ 1050C 1051 INTEGER NRAND, MAX 1052 DOUBLE PRECISION FA01AS 1053 EXTERNAL FA01AS 1054 NRAND=INT(FA01AS(1)*FLOAT(MAX))+1 1055 RETURN 1056 END 1057C 1058 SUBROUTINE FA01CS(IL,IR) 1059C ======================== 1060C 1061 DOUBLE PRECISION G 1062 INTEGER IL,IR 1063 COMMON/FA01ES/G 1064 G= 1431655765.D0 1065C 1066 IL=G/65536.D0 1067 IR=G-65536.D0*FLOAT(IL) 1068 RETURN 1069 END 1070C 1071 SUBROUTINE FA01DS(IL,IR) 1072C ======================== 1073C 1074 DOUBLE PRECISION G 1075 INTEGER IL,IR 1076 COMMON/FA01ES/G 1077 G= 1431655765.D0 1078C 1079 G=65536.D0*FLOAT(IL)+FLOAT(IR) 1080 RETURN 1081 END 1082C 1083C_BEGIN_FM02AD 1084C 1085 DOUBLE PRECISION FUNCTION FM02AD(N,A,IA,B,IB) 1086C ============================================= 1087C 1088C Compute the inner product of two double precision real 1089C vectors accumulating the result double precision, when the 1090C elements of each vector are stored at some fixed displacement 1091C from neighbouring elements. Given vectors A={a(j)}, 1092C B={b(j)} of length N, evaluates w=a(j)b(j) summed over 1093C j=1..N. Can be used to evaluate inner products involving 1094C rows of multi-dimensional arrays. 1095C It can be used as an alternative to the assembler version, 1096C but note that it is likely to be significantly slower in execution. 1097C 1098C Argument list 1099C ------------- 1100C 1101C N (I) (integer) The length of the vectors (if N <= 0 FM02AD = 0) 1102C 1103C A (I) (double precision) The first vector 1104C 1105C IA (I) (integer) Subscript displacement between elements of A 1106C 1107C B (I) (double precision) The second vector 1108C 1109C IB (I) (integer) Subscript displacement between elements of B 1110C 1111C FM02AD the result 1112C 1113C 1114C_END_FM02AD 1115C 1116 DOUBLE PRECISION R1,A,B 1117C 1118C---- The following statement changed from a(n),b(n) to avoid vax dynamic 1119C array check failure. 1120C 1121 DIMENSION A(1),B(1) 1122 INTEGER N,JA,IA,JB,IB,I 1123C 1124 R1=0D0 1125 IF(N.LE.0) GO TO 2 1126 JA=1 1127 IF(IA.LT.0) JA=1-(N-1)*IA 1128 JB=1 1129 IF(IB.LT.0) JB=1-(N-1)*IB 1130 I=0 1131 1 I=I+1 1132 R1=R1+A(JA)*B(JB) 1133 JA=JA+IA 1134 JB=JB+IB 1135 IF(I.LT.N) GO TO 1 1136 2 FM02AD=R1 1137 RETURN 1138C 1139 END 1140C 1141C 1142C_BEGIN_ICROSS 1143C 1144 SUBROUTINE ICROSS(A,B,C) 1145C ======================== 1146C 1147C Cross product (integer version) 1148C 1149C A = B x C 1150C 1151C .. Array Arguments .. 1152 INTEGER A(3),B(3),C(3) 1153C 1154C_END_ICROSS 1155C .. 1156 A(1) = B(2)*C(3) - C(2)*B(3) 1157 A(2) = B(3)*C(1) - C(3)*B(1) 1158 A(3) = B(1)*C(2) - C(1)*B(2) 1159 END 1160C 1161C 1162C_BEGIN_IDOT 1163C 1164 INTEGER FUNCTION IDOT(A,B) 1165C ========================== 1166C 1167C Dot product (integer version) 1168C 1169C IDOT = A . B 1170C 1171C .. Array Arguments .. 1172 INTEGER A(3),B(3) 1173C 1174C_END_IDOT 1175C .. 1176 IDOT = A(1)*B(1) + A(2)*B(2) + A(3)*B(3) 1177 END 1178C 1179C 1180C_BEGIN_IMINV3 1181C 1182 SUBROUTINE IMINV3(A,B,D) 1183C ======================= 1184C 1185C Invert a general 3x3 matrix and return determinant in D 1186C (integer version) 1187C 1188C A = (B)-1 1189C 1190C .. Scalar Arguments .. 1191 INTEGER D 1192C .. 1193C .. Array Arguments .. 1194 INTEGER A(3,3),B(3,3) 1195C 1196C_END_IMINV3 1197C .. 1198C .. Local Scalars .. 1199 INTEGER I,J 1200C .. 1201C .. Local Arrays .. 1202 INTEGER C(3,3) 1203C .. 1204C .. External Functions .. 1205 INTEGER IDOT 1206 EXTERNAL IDOT 1207C .. 1208C .. External Subroutines .. 1209 EXTERNAL ICROSS 1210C .. 1211C .. Intrinsic Functions .. 1212 INTRINSIC ABS 1213C .. 1214 CALL ICROSS(C(1,1),B(1,2),B(1,3)) 1215 CALL ICROSS(C(1,2),B(1,3),B(1,1)) 1216 CALL ICROSS(C(1,3),B(1,1),B(1,2)) 1217 D = IDOT(B(1,1),C(1,1)) 1218C 1219C---- Test determinant 1220C 1221 IF (ABS(D).GT.0) THEN 1222C 1223C---- Determinant is non-zero 1224C 1225 DO 20 I = 1,3 1226 DO 10 J = 1,3 1227 A(I,J) = C(J,I)/D 1228 10 CONTINUE 1229 20 CONTINUE 1230 1231 ELSE 1232 D = 0 1233 END IF 1234 1235 END 1236C 1237C 1238C_BEGIN_MATMUL 1239C 1240 SUBROUTINE MATMUL(A,B,C) 1241C ======================== 1242C 1243C Multiply two 3x3 matrices 1244C 1245C A = BC 1246C 1247C .. Array Arguments .. 1248 REAL A(3,3),B(3,3),C(3,3) 1249C 1250C_END_MATMUL 1251C .. 1252C .. Local Scalars .. 1253 INTEGER I,J,K 1254C .. 1255 DO 30 I = 1,3 1256 DO 20 J = 1,3 1257 A(I,J) = 0.0 1258 DO 10 K = 1,3 1259 A(I,J) = B(I,K)*C(K,J) + A(I,J) 1260 10 CONTINUE 1261 20 CONTINUE 1262 30 CONTINUE 1263 END 1264C 1265C 1266C_BEGIN_MATMULNM 1267C 1268 SUBROUTINE MATMULNM(N,M,A,B,C) 1269C ======================== 1270C 1271C Multiply NxM MXN matrices 1272C 1273C A = BC 1274C 1275C .. Array Arguments .. 1276 INTEGER N,M 1277 REAL A(N,N),B(N,M),C(M,N) 1278C 1279C_END_MATMUL 1280C .. 1281C .. Local Scalars .. 1282 INTEGER I,J,K 1283C .. 1284 DO 30 I = 1,N 1285 DO 20 J = 1,N 1286 A(I,J) = 0.0 1287 DO 10 K = 1,M 1288 A(I,J) = B(I,K)*C(K,J) + A(I,J) 1289 10 CONTINUE 1290 20 CONTINUE 1291 30 CONTINUE 1292 END 1293C 1294C 1295C_BEGIN_MATVEC 1296C 1297 SUBROUTINE MATVEC(V,A,B) 1298C ======================== 1299C 1300C Post-multiply a 3x3 matrix by a vector 1301C 1302C V = AB 1303C 1304C .. Array Arguments .. 1305 REAL A(3,3),B(3),V(3) 1306C 1307C_END_MATVEC 1308C .. 1309C .. Local Scalars .. 1310 REAL S 1311 INTEGER I,J 1312C .. 1313 DO 20 I = 1,3 1314 S = 0 1315 DO 10 J = 1,3 1316 S = A(I,J)*B(J) + S 1317 10 CONTINUE 1318 V(I) = S 1319 20 CONTINUE 1320 END 1321C 1322C 1323C_BEGIN_MATMULGEN 1324C 1325 SUBROUTINE MATMULGEN(Nb,Mbc,Nc,A,B,C) 1326C ===================================== 1327C 1328C Generalised matrix multiplication subroutine 1329C Multiplies a NbxMbc matrix (B) by a MbcXNc 1330C (C) matrix, so that 1331C 1332C A = BC 1333C 1334 IMPLICIT NONE 1335C .. 1336C .. Scalar Arguments .. 1337 INTEGER Nb,Mbc,Nc 1338C .. 1339C .. Array Arguments .. 1340 REAL A(Nb,Nc),B(Nb,Mbc),C(Mbc,Nc) 1341C 1342C_END_MATMULGEN 1343C .. 1344C .. Local Scalars .. 1345 INTEGER I,J,K 1346C .. 1347 DO 30 J = 1,Nc 1348 DO 20 I = 1,Nb 1349 A(I,J) = 0.0 1350 DO 10 K = 1,Mbc 1351 A(I,J) = B(I,K)*C(K,J) + A(I,J) 1352 10 CONTINUE 1353 20 CONTINUE 1354 30 CONTINUE 1355 END 1356C 1357C 1358C_BEGIN_MC04B 1359C 1360 SUBROUTINE MC04B(A,ALPHA,BETA,M,IA,Q) 1361C ===================================== 1362C 1363C Transforms a real symmetric matrix A={a(i,j)}, i, j=1..IA 1364C into a tri-diagonal matrix having the same eigenvalues as A 1365C using Householder's method. 1366C 1367C .. Scalar Arguments .. 1368 INTEGER IA,M 1369C .. 1370C .. Array Arguments .. 1371 REAL A(IA,1),ALPHA(1),BETA(1),Q(1) 1372C 1373C_END_MC04B 1374C .. 1375C .. Local Scalars .. 1376 REAL BIGK,H,PP,PP1,QJ 1377 INTEGER I,I1,I2,J,J1,KI,KJ,M1,M2 1378C .. 1379C .. Intrinsic Functions .. 1380 INTRINSIC SQRT 1381C .. 1382 ALPHA(1) = A(1,1) 1383 DO 20 J = 2,M 1384 J1 = J - 1 1385 DO 10 I = 1,J1 1386 A(I,J) = A(J,I) 1387 10 CONTINUE 1388 ALPHA(J) = A(J,J) 1389 20 CONTINUE 1390 M1 = M - 1 1391 M2 = M - 2 1392 DO 110 I = 1,M2 1393 PP = 0.0 1394 I1 = I + 1 1395 DO 30 J = I1,M 1396 PP = A(I,J)**2 + PP 1397 30 CONTINUE 1398 PP1 = SQRT(PP) 1399 IF (A(I,I+1).LT.0) THEN 1400 BETA(I+1) = PP1 1401 1402 ELSE 1403 BETA(I+1) = -PP1 1404 END IF 1405 1406 IF (PP.GT.0) THEN 1407 H = PP - BETA(I+1)*A(I,I+1) 1408 A(I,I+1) = A(I,I+1) - BETA(I+1) 1409 DO 60 KI = I1,M 1410 QJ = 0.0 1411 DO 40 KJ = I1,KI 1412 QJ = A(KJ,KI)*A(I,KJ) + QJ 1413 40 CONTINUE 1414 IF (KI-M.LT.0) THEN 1415 I2 = KI + 1 1416 DO 50 KJ = I2,M 1417 QJ = A(KI,KJ)*A(I,KJ) + QJ 1418 50 CONTINUE 1419 END IF 1420 1421 Q(KI) = QJ/H 1422 60 CONTINUE 1423 BIGK = 0.0 1424 DO 70 KJ = I1,M 1425 BIGK = A(I,KJ)*Q(KJ) + BIGK 1426 70 CONTINUE 1427 BIGK = BIGK/ (2.0*H) 1428 DO 80 KJ = I1,M 1429 Q(KJ) = Q(KJ) - A(I,KJ)*BIGK 1430 80 CONTINUE 1431 DO 100 KI = I1,M 1432 DO 90 KJ = KI,M 1433 A(KI,KJ) = A(KI,KJ) - Q(KI)*A(I,KJ) - 1434 + Q(KJ)*A(I,KI) 1435 90 CONTINUE 1436 100 CONTINUE 1437 END IF 1438 1439 110 CONTINUE 1440 DO 120 I = 2,M 1441 H = ALPHA(I) 1442 ALPHA(I) = A(I,I) 1443 A(I,I) = H 1444 120 CONTINUE 1445 BETA(M) = A(M-1,M) 1446 END 1447C 1448C 1449C_BEGIN_MINVN 1450C 1451 SUBROUTINE MINVN(A,N,D,L,M) 1452C =========================== 1453C 1454C 1455C---- Purpose 1456C ======= 1457C 1458C invert a matrix 1459C 1460C---- Usage 1461C ====== 1462C 1463C CALL MINVN(A,N,D,L,M) 1464C 1465C---- Description of parameters 1466C ========================= 1467C 1468C A - input matrix, destroyed in computation and replaced by 1469C resultant inverse. 1470C 1471C N - order of matrix A 1472C 1473C D - resultant determinant 1474C 1475C L - work vector of length n 1476C 1477C M - work vector of length n 1478C 1479C---- Remarks 1480C ======= 1481C 1482C Matrix a must be a general matrix 1483C 1484C---- Subroutines and function subprograms required 1485C ============================================= 1486C 1487C NONE 1488C 1489C---- Method 1490C ====== 1491C 1492C The standard gauss-jordan method is used. the determinant 1493C is also calculated. a determinant of zero indicates that 1494C the matrix is singular. 1495C 1496C 1497C---- Note 1498C ===== 1499C 1500C If a double precision version of this routine is desired, the 1501C c in column 1 should be removed from the double precision 1502C statement which follows. 1503C 1504C double precision a,d,biga,hold 1505C 1506C the c must also be removed from double precision statements 1507C appearing in other routines used in conjunction with this 1508C routine. 1509C 1510C The double precision version of this subroutine must also 1511C contain double precision fortran functions. abs in statement 1512C 10 must be changed to dabs. 1513C 1514ccc REAL*8 D 1515C 1516C_END_MINVN 1517C 1518C---- Search for largest element 1519C 1520C .. Scalar Arguments .. 1521 REAL D 1522 INTEGER N 1523C .. 1524C .. Array Arguments .. 1525 REAL A(N*N) 1526 INTEGER L(N),M(N) 1527C .. 1528C .. Local Scalars .. 1529 REAL BIGA,HOLD 1530 INTEGER I,IJ,IK,IZ,J,JI,JK,JP,JQ,JR,K,KI,KJ,KK,NK 1531C .. 1532C .. Intrinsic Functions .. 1533 INTRINSIC ABS 1534C .. 1535C 1536C 1537 D = 1.0 1538 NK = -N 1539 DO 90 K = 1,N 1540 NK = NK + N 1541 L(K) = K 1542 M(K) = K 1543 KK = NK + K 1544 BIGA = A(KK) 1545 DO 20 J = K,N 1546 IZ = (J-1)*N 1547 DO 10 I = K,N 1548 IJ = IZ + I 1549 IF ((ABS(BIGA)-ABS(A(IJ))).LT.0.0) THEN 1550 BIGA = A(IJ) 1551 L(K) = I 1552 M(K) = J 1553 END IF 1554 10 CONTINUE 1555 20 CONTINUE 1556C 1557C---- Interchange rows 1558C 1559 J = L(K) 1560 IF ((J-K).GT.0) THEN 1561 KI = K - N 1562 DO 30 I = 1,N 1563 KI = KI + N 1564 HOLD = -A(KI) 1565 JI = KI - K + J 1566 A(KI) = A(JI) 1567 A(JI) = HOLD 1568 30 CONTINUE 1569 END IF 1570C 1571C---- Interchange columns 1572C 1573 I = M(K) 1574 IF ((I-K).GT.0) THEN 1575 JP = (I-1)*N 1576 DO 40 J = 1,N 1577 JK = NK + J 1578 JI = JP + J 1579 HOLD = -A(JK) 1580 A(JK) = A(JI) 1581 A(JI) = HOLD 1582 40 CONTINUE 1583 END IF 1584C 1585C---- Divide column by minus pivot (value of pivot element is 1586C contained in biga) 1587C 1588 IF (BIGA.NE.0.0) THEN 1589 DO 50 I = 1,N 1590 IF ((I-K).NE.0) THEN 1591 IK = NK + I 1592 A(IK) = A(IK)/ (-BIGA) 1593 END IF 1594 50 CONTINUE 1595C 1596C---- Reduce matrix 1597C 1598 DO 70 I = 1,N 1599 IK = NK + I 1600 HOLD = A(IK) 1601 IJ = I - N 1602 DO 60 J = 1,N 1603 IJ = IJ + N 1604 IF ((I-K).NE.0) THEN 1605 IF ((J-K).NE.0) THEN 1606 KJ = IJ - I + K 1607 A(IJ) = A(KJ)*HOLD + A(IJ) 1608 END IF 1609 END IF 1610 60 CONTINUE 1611 70 CONTINUE 1612C 1613C---- Divide row by pivot 1614C 1615 KJ = K - N 1616 DO 80 J = 1,N 1617 KJ = KJ + N 1618 IF ((J-K).NE.0) A(KJ) = A(KJ)/BIGA 1619 80 CONTINUE 1620C 1621C---- Product of pivots 1622C 1623 D = D*BIGA 1624C 1625C---- Replace pivot by reciprocal 1626C 1627 A(KK) = 1.0/BIGA 1628 ELSE 1629 GO TO 130 1630 END IF 1631 90 CONTINUE 1632C 1633C---- Final row and column interchange 1634C 1635 K = N 1636 100 CONTINUE 1637 K = (K-1) 1638 IF (K.GT.0) THEN 1639 I = L(K) 1640 IF ((I-K).GT.0) THEN 1641 JQ = (K-1)*N 1642 JR = (I-1)*N 1643 DO 110 J = 1,N 1644 JK = JQ + J 1645 HOLD = A(JK) 1646 JI = JR + J 1647 A(JK) = -A(JI) 1648 A(JI) = HOLD 1649 110 CONTINUE 1650 END IF 1651 J = M(K) 1652 IF ((J-K).GT.0) THEN 1653 KI = K - N 1654 DO 120 I = 1,N 1655 KI = KI + N 1656 HOLD = A(KI) 1657 JI = KI - K + J 1658 A(KI) = -A(JI) 1659 A(JI) = HOLD 1660 120 CONTINUE 1661 END IF 1662 GO TO 100 1663 ELSE 1664 RETURN 1665 END IF 1666 130 D = 0.0 1667C 1668C 1669 END 1670C 1671C 1672C_BEGIN_MINV3 1673C 1674 SUBROUTINE MINV3(A,B,D) 1675C ====================== 1676C 1677C Invert a general 3x3 matrix and return determinant in D 1678C 1679C A = (B)-1 1680C 1681C .. Scalar Arguments .. 1682 REAL D 1683C .. 1684C .. Array Arguments .. 1685 REAL A(3,3),B(3,3) 1686C 1687C_END_MINV3 1688C .. 1689C .. Local Scalars .. 1690 INTEGER I,J 1691C .. 1692C .. Local Arrays .. 1693 REAL C(3,3) 1694C .. 1695C .. External Functions .. 1696 REAL DOT 1697 EXTERNAL DOT 1698C .. 1699C .. External Subroutines .. 1700 EXTERNAL CROSS 1701C .. 1702C .. Intrinsic Functions .. 1703 INTRINSIC ABS 1704C .. 1705 CALL CROSS(C(1,1),B(1,2),B(1,3)) 1706 CALL CROSS(C(1,2),B(1,3),B(1,1)) 1707 CALL CROSS(C(1,3),B(1,1),B(1,2)) 1708 D = DOT(B(1,1),C(1,1)) 1709C 1710C---- Test determinant 1711C 1712 IF (ABS(D).GT.1.0E-30) THEN 1713C 1714C---- Determinant is non-zero 1715C 1716 DO 20 I = 1,3 1717 DO 10 J = 1,3 1718 A(I,J) = C(J,I)/D 1719 10 CONTINUE 1720 20 CONTINUE 1721 1722 ELSE 1723 D = 0.0 1724 END IF 1725 1726 END 1727C 1728C 1729C_BEGIN_RANMAR 1730C 1731 SUBROUTINE RANMAR(RVEC,LEN) 1732C =========================== 1733C 1734C Universal random number generator proposed by Marsaglia and Zaman 1735C in report FSU-SCRI-87-50 1736C slightly modified by F. James, 1988 to generate a vector 1737C of pseudorandom numbers RVEC of length LEN 1738C and making the COMMON block include everything needed to 1739C specify completely the state of the generator. 1740C Transcribed from CERN report DD/88/22. 1741C Rather inelegant messing about added by D. Love, Jan. 1989 to 1742C make sure initialisation always occurs. 1743C *** James says that this is the preferred generator. 1744C Gives bit-identical results on all machines with at least 1745C 24-bit mantissas in the flotaing point representation (i.e. 1746C all common 32-bit computers. Fairly fast, satisfies very 1747C stringent tests, has very long period and makes it very 1748C simple to generate independly disjoint sequences. 1749C See also RANECU. 1750C The state of the generator may be saved/restored using the 1751C whole contents of /RASET1/. 1752C Call RANMAR to get a vector, RMARIN to initialise. 1753C 1754C Argument list 1755C ------------- 1756C 1757C VREC (O) (REAL) Random Vector 1758C 1759C LEN (I) (INTEGER) Length of random vector 1760C 1761C 1762C For ENTRY point RMARIN 1763C ---------------------- 1764C 1765C Initialisation for RANMAR. The input values should 1766C be in the ranges: 0<=ij<=31328, 0<=kl<=30081 1767C This shows the correspondence between the simplified input seeds 1768C IJ, KL and the original Marsaglia-Zaman seeds i,j,k,l 1769C To get standard values in Marsaglia-Zaman paper, 1770C (I=12, J=34, K=56, L=78) put IJ=1802, KL=9373 1771C 1772C IJ (I) (INTEGER) Seed for random number generator 1773C 1774C KL (I) (INTEGER) Seed for randon number generator 1775C 1776C_END_RANMAR 1777C 1778C .. 1779C .. Agruments .. 1780 REAL RVEC(*) 1781 INTEGER LEN,IJ,KL 1782C .. 1783C .. Common Variables .. 1784 REAL C,CD,CM,U 1785 INTEGER I97,J97 1786C .. 1787C .. Local Scalars .. 1788 REAL S,T,UNI 1789 INTEGER I,II,IVEC,J,JJ,K,L,M 1790 LOGICAL INITED 1791C .. 1792C .. Intrinsic Functions .. 1793 INTRINSIC MOD 1794C .. 1795C .. Common Blocks .. 1796 COMMON /RASET1/ U(97),C,CD,CM,I97,J97 1797C .. 1798C .. Save Statement .. 1799 SAVE INITED, /RASET1/ 1800C .. 1801C .. Data Statement .. 1802 DATA INITED /.FALSE./ 1803C 1804C---- If initialised, fill RVEC and RETURN. If not, do initialisation 1805C and return here later. 1806C 1807 1 IF (INITED) THEN 1808 DO 100 IVEC=1,LEN 1809 UNI=U(I97)-U(J97) 1810 IF (UNI.LT.0.) UNI=UNI+1. 1811 U(I97)=UNI 1812 I97=I97-1 1813 IF (I97.EQ.0) I97=97 1814 J97=J97-1 1815 IF (J97.EQ.0) J97=97 1816 C=C-CD 1817 IF (C.LT.0.) C=C+CM 1818 UNI=UNI-C 1819 IF (UNI.LT.0.) UNI=UNI+1. 1820 RVEC(IVEC)=UNI 1821 100 CONTINUE 1822 RETURN 1823 ENDIF 1824 I=MOD(1802/177,177)+2 1825 J=MOD(1802,177)+2 1826 K=MOD(9373/169,178)+1 1827 L=MOD(9373,169) 1828C 1829 GOTO 10 1830C 1831C---- Initialise and return without filling RVEC 1832C 1833 ENTRY RMARIN(IJ,KL) 1834 I=MOD(IJ/177,177)+2 1835 J=MOD(IJ,177)+2 1836 K=MOD(KL/169,178)+1 1837 L=MOD(KL,169) 1838 INITED=.TRUE. 1839 1840 10 CONTINUE 1841 DO 2 II=1,97 1842 S=0. 1843 T=.5 1844 DO 3 JJ=1,24 1845 M=MOD(MOD(I*J,179)*K,179) 1846 I=J 1847 J=K 1848 K=M 1849 L=MOD(53*L+1,169) 1850 IF (MOD(L*M,64).GE.32) S=S+T 1851 T=0.5*T 1852 3 CONTINUE 1853 U(II)=S 1854 2 CONTINUE 1855 C=362436./16777216. 1856 CD=7654321./16777216. 1857 CM=16777213./16777216. 1858 I97=97 1859 J97=33 1860 IF (.NOT. INITED) THEN 1861 INITED=.TRUE. 1862 GOTO 1 1863 ENDIF 1864 1865 END 1866C 1867C 1868C_BEGIN_SCALEV 1869C 1870 SUBROUTINE SCALEV(A,X,B) 1871C ======================== 1872C 1873C Scale vector B with scalar X and put result in A 1874C 1875C .. Scalar Arguments .. 1876 REAL X 1877C .. 1878C .. Array Arguments .. 1879 REAL A(3),B(3) 1880C 1881C_END_SCALEV 1882C .. 1883C .. Local Scalars .. 1884 INTEGER I 1885C .. 1886 DO 10 I = 1,3 1887 A(I) = B(I)*X 1888 10 CONTINUE 1889 END 1890C 1891C 1892C_BEGIN_TRANSP 1893C 1894 SUBROUTINE TRANSP(A,B) 1895C ====================== 1896C 1897C---- Transpose a 3x3 matrix 1898C 1899C A = BT 1900C 1901C .. Array Arguments .. 1902 REAL A(3,3),B(3,3) 1903C 1904C_END_TRANSP 1905C .. 1906C .. Local Scalars .. 1907 INTEGER I,J 1908C .. 1909 DO 20 I = 1,3 1910 DO 10 J = 1,3 1911 A(I,J) = B(J,I) 1912 10 CONTINUE 1913 20 CONTINUE 1914 END 1915C 1916C 1917C_BEGIN_UNIT 1918C 1919 SUBROUTINE UNIT(V) 1920C ================= 1921C 1922C Vector V reduced to unit vector 1923C 1924C .. Array Arguments .. 1925 REAL V(3) 1926C 1927C_END_UNIT 1928C .. 1929C .. Local Scalars .. 1930 REAL VMOD 1931 INTEGER I 1932C .. 1933C .. Intrinsic Functions .. 1934 INTRINSIC SQRT 1935C .. 1936 VMOD = V(1)**2 + V(2)**2 + V(3)**2 1937 VMOD = SQRT(VMOD) 1938 DO 10 I = 1,3 1939 V(I) = V(I)/VMOD 1940 10 CONTINUE 1941 END 1942C 1943C 1944C_BEGIN_VDIF 1945C 1946 SUBROUTINE VDIF(A,B,C) 1947C ===================== 1948C 1949C Subtract two vectors 1950C 1951C A = B - C 1952C 1953C .. Array Arguments .. 1954 REAL A(3),B(3),C(3) 1955C 1956C_END_VDIF 1957C .. 1958C .. Local Scalars .. 1959 INTEGER I 1960C .. 1961 DO 10 I = 1,3 1962 A(I) = B(I) - C(I) 1963 10 CONTINUE 1964 END 1965C 1966C 1967C_BEGIN_VSET 1968C 1969 SUBROUTINE VSET(A,B) 1970C ==================== 1971C 1972C Copy a vector from B to A 1973C 1974C .. Array Arguments .. 1975 REAL A(3),B(3) 1976C 1977C_END_VSET 1978C .. 1979C .. Local Scalars .. 1980 INTEGER I 1981C .. 1982 DO 10 I = 1,3 1983 A(I) = B(I) 1984 10 CONTINUE 1985 END 1986C 1987C 1988C_BEGIN_VSUM 1989C 1990 SUBROUTINE VSUM(A,B,C) 1991C ====================== 1992C 1993C Add two vectors 1994C 1995C A = B + C 1996C 1997C .. Array Arguments .. 1998 REAL A(3),B(3),C(3) 1999C 2000C_END_VSUM 2001C .. 2002C .. Local Scalars .. 2003 INTEGER I 2004C .. 2005 DO 10 I = 1,3 2006 A(I) = B(I) + C(I) 2007 10 CONTINUE 2008 END 2009C 2010C 2011C_BEGIN_ZIPIN 2012C 2013 SUBROUTINE ZIPIN(ID,N,BUF) 2014C ========================== 2015C 2016C Fast binary read on unit ID into real array BUF of length N 2017C 2018C .. Scalar Arguments .. 2019 INTEGER ID,N 2020C .. 2021C .. Array Arguments .. 2022 REAL BUF(N) 2023C 2024C_END_ZIPIN 2025C .. 2026 READ (ID) BUF 2027 END 2028C 2029C 2030C_BEGIN_ZIPOUT 2031C 2032 SUBROUTINE ZIPOUT(ID,N,BUF) 2033C =========================== 2034C 2035C Fast binary write to unit ID of real array BUF length N 2036C 2037C .. Scalar Arguments .. 2038 INTEGER ID,N 2039C .. 2040C .. Array Arguments .. 2041 REAL BUF(N) 2042C 2043C_END_ZIPOUT 2044C .. 2045 WRITE (ID) BUF 2046 END 2047C 2048C The following routines calculate Bessel functions of 2049C the first kind, their derivatives, and their zeros. 2050C Subroutine ZJVX is from Ian Tickle. 2051C Subroutines JDVX, JVX, GAMMA and functions MSTA1, MSTA2 2052C are from MJYV at http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html 2053C These routines are copyrighted, but permission is given to incorporate these 2054C routines into programs provided that the copyright is acknowledged. 2055 2056 SUBROUTINE ZJVX(V,IZ,BJ,DJ,X) 2057C 2058C ======================================================= 2059C Purpose: Compute zero of Bessel function Jv(x) 2060C Input : V --- Order of Jv(x) 2061C IZ --- Index of previous zero of Jv(x) 2062C X --- Value of previous zero of Jv(x) 2063C Output: BJ(N) --- Jv(x) for N = 0...INT(V) 2064C DJ(N) --- J'v(x) for N = 0...INT(V) 2065C X --- Value of next zero of Jv(x) 2066C 2067 IMPLICIT NONE 2068 INTEGER IZ,N 2069 REAL*8 V,VM,X,X0 2070 REAL*8 BJ(0:*),DJ(0:*) 2071C 2072C#### Initial guess at zero: needs previous value if not first. 2073 N=V 2074 IF (IZ.EQ.0) THEN 2075 X=1.99535+.8333883*SQRT(V)+.984584*V 2076 ELSEIF (N.LE.10) THEN 2077 X=X+3.11+.0138*V+(.04832+.2804*V)/(IZ+1)**2 2078 ELSE 2079 X=X+3.001+.0105*V+(11.52+.48525*V)/(IZ+3)**2 2080 ENDIF 2081C 2082C#### Polish zero by Newton-Cotes iteration. 20831 CALL JDVX(V,X,VM,BJ,DJ) 2084 IF (INT(VM).NE.N) CALL CCPERR(1,'VM != N in ZJVX.') 2085 X0=X 2086 X=X-BJ(N)/DJ(N) 2087 IF (ABS(X-X0).GT.1D-10) GOTO 1 2088 END 2089C 2090C 2091 SUBROUTINE JDVX(V,X,VM,BJ,DJ) 2092C 2093C ======================================================= 2094C Purpose: Compute Bessel functions Jv(x) and their derivatives 2095C Input : x --- Argument of Jv(x) 2096C v --- Order of Jv(x) 2097C ( v = n+v0, 0 <= v0 < 1, n = 0,1,2,... ) 2098C Output: BJ(n) --- Jn+v0(x) 2099C DJ(n) --- Jn+v0'(x) 2100C VM --- Highest order computed 2101C Routines called: 2102C (1) GAMMA for computing gamma function 2103C (2) MSTA1 and MSTA2 for computing the starting 2104C point for backward recurrence 2105C ======================================================= 2106C 2107 IMPLICIT NONE 2108 INTEGER J,K,K0,L,M,N 2109 REAL*8 A,A0,BJV0,BJV1,BJVL,CK,CS,F,F0,F1,F2,GA,PI,PX,QX,R,RP,RP2, 2110 &RQ,SK,V,V0,VG,VL,VM,VV,X,X2,XK 2111 PARAMETER(PI=3.141592653589793D0, RP2=.63661977236758D0) 2112 REAL*8 BJ(0:*),DJ(0:*) 2113 INTEGER MSTA1,MSTA2 2114C 2115 X2=X*X 2116 N=INT(V) 2117 V0=V-N 2118 IF (X.LT.1D-100) THEN 2119 DO K=0,N 2120 BJ(K)=0D0 2121 DJ(K)=0D0 2122 ENDDO 2123 IF (V0.EQ.0D0) THEN 2124 BJ(0)=1D0 2125 DJ(1)=.5D0 2126 ELSE 2127 DJ(0)=1D300 2128 ENDIF 2129 VM=V 2130 ELSE 2131 IF (X.LE.12D0) THEN 2132 DO L=0,1 2133 VL=V0+L 2134 BJVL=1D0 2135 R=1D0 2136 DO K=1,40 2137 R=-.25D0*R*X2/(K*(K+VL)) 2138 BJVL=BJVL+R 2139 IF (ABS(R).LT.ABS(BJVL)*1D-15) GOTO 20 2140 ENDDO 214120 VG=1D0+VL 2142 CALL GAMMA(VG,GA) 2143 A=(.5D0*X)**VL/GA 2144 IF (L.EQ.0) THEN 2145 BJV0=BJVL*A 2146 ELSEIF (L.EQ.1) THEN 2147 BJV1=BJVL*A 2148 ENDIF 2149 ENDDO 2150 ELSE 2151 K0=11 2152 IF (X.GE.35D0) K0=10 2153 IF (X.GE.50D0) K0=8 2154 DO J=0,1 2155 VV=4D0*(J+V0)*(J+V0) 2156 PX=1D0 2157 RP=1D0 2158 DO K=1,K0 2159 RP=-.78125D-2*RP*(VV-(4*K-3)**2)*(VV-(4*K-1)**2)/ 2160 & (K*(2*K-1)*X2) 2161 PX=PX+RP 2162 ENDDO 2163 QX=1D0 2164 RQ=1D0 2165 DO K=1,K0 2166 RQ=-.78125D-2*RQ*(VV-(4*K-1)**2)*(VV-(4*K+1)**2)/ 2167 & (K*(2*K+1)*X2) 2168 QX=QX+RQ 2169 ENDDO 2170 QX=.125D0*(VV-1D0)*QX/X 2171 XK=X-(.5D0*(J+V0)+.25D0)*PI 2172 A0=SQRT(RP2/X) 2173 CK=COS(XK) 2174 SK=SIN(XK) 2175 IF (J.EQ.0) THEN 2176 BJV0=A0*(PX*CK-QX*SK) 2177 ELSEIF (J.EQ.1) THEN 2178 BJV1=A0*(PX*CK-QX*SK) 2179 ENDIF 2180 ENDDO 2181 ENDIF 2182 BJ(0)=BJV0 2183 BJ(1)=BJV1 2184 DJ(0)=V0/X*BJ(0)-BJ(1) 2185 DJ(1)=BJ(0)-(1D0+V0)/X*BJ(1) 2186 IF (N.GE.2 .AND. N.LE.INT(.9*X)) THEN 2187 F0=BJV0 2188 F1=BJV1 2189 DO K=2,N 2190 F=2D0*(K+V0-1D0)/X*F1-F0 2191 BJ(K)=F 2192 F0=F1 2193 F1=F 2194 ENDDO 2195 ELSEIF (N.GE.2) THEN 2196 M=MSTA1(X,200) 2197 IF (M.LT.N) THEN 2198 N=M 2199 ELSE 2200 M=MSTA2(X,N,15) 2201 ENDIF 2202 F2=0D0 2203 F1=1D-100 2204 DO K=M,0,-1 2205 F=2D0*(V0+K+1D0)/X*F1-F2 2206 IF (K.LE.N) BJ(K)=F 2207 F2=F1 2208 F1=F 2209 ENDDO 2210 IF (ABS(BJV0).GT.ABS(BJV1)) THEN 2211 CS=BJV0/F 2212 ELSE 2213 CS=BJV1/F2 2214 ENDIF 2215 DO K=0,N 2216 BJ(K)=CS*BJ(K) 2217 ENDDO 2218 ENDIF 2219 DO K=2,N 2220 DJ(K)=BJ(K-1)-(K+V0)/X*BJ(K) 2221 ENDDO 2222 VM=N+V0 2223 ENDIF 2224 END 2225C 2226C 2227 SUBROUTINE JVX(V,X,VM,BJ) 2228C 2229C ======================================================= 2230C Purpose: Compute Bessel functions Jv(x) 2231C Input : x --- Argument of Jv(x) 2232C v --- Order of Jv(x) 2233C ( v = n+v0, 0 <= v0 < 1, n = 0,1,2,... ) 2234C Output: BJ(n) --- Jn+v0(x) 2235C VM --- Highest order computed 2236C Routines called: 2237C (1) GAMMA for computing gamma function 2238C (2) MSTA1 and MSTA2 for computing the starting 2239C point for backward recurrence 2240C ======================================================= 2241C 2242 IMPLICIT NONE 2243 INTEGER J,K,K0,L,M,N 2244 REAL*8 A,A0,BJV0,BJV1,BJVL,CK,CS,F,F0,F1,F2,GA,PI,PX,QX,R,RP,RP2, 2245 &RQ,SK,V,V0,VG,VL,VM,VV,X,X2,XK 2246 PARAMETER(PI=3.141592653589793D0, RP2=.63661977236758D0) 2247 REAL*8 BJ(0:*) 2248 INTEGER MSTA1,MSTA2 2249C 2250 X2=X*X 2251 N=INT(V) 2252 V0=V-N 2253 IF (X.LT.1D-100) THEN 2254 DO K=0,N 2255 BJ(K)=0D0 2256 ENDDO 2257 IF (V0.EQ.0D0) BJ(0)=1D0 2258 VM=V 2259 ELSE 2260 IF (X.LE.12D0) THEN 2261 DO L=0,1 2262 VL=V0+L 2263 BJVL=1D0 2264 R=1D0 2265 DO K=1,40 2266 R=-.25D0*R*X2/(K*(K+VL)) 2267 BJVL=BJVL+R 2268 IF (ABS(R).LT.ABS(BJVL)*1D-15) GOTO 20 2269 ENDDO 227020 VG=1D0+VL 2271 CALL GAMMA(VG,GA) 2272 A=(.5D0*X)**VL/GA 2273 IF (L.EQ.0) THEN 2274 BJV0=BJVL*A 2275 ELSEIF (L.EQ.1) THEN 2276 BJV1=BJVL*A 2277 ENDIF 2278 ENDDO 2279 ELSE 2280 K0=11 2281 IF (X.GE.35D0) K0=10 2282 IF (X.GE.50D0) K0=8 2283 DO J=0,1 2284 VV=4D0*(J+V0)*(J+V0) 2285 PX=1D0 2286 RP=1D0 2287 DO K=1,K0 2288 RP=-.78125D-2*RP*(VV-(4*K-3)**2)*(VV-(4*K-1)**2)/ 2289 & (K*(2*K-1)*X2) 2290 PX=PX+RP 2291 ENDDO 2292 QX=1D0 2293 RQ=1D0 2294 DO K=1,K0 2295 RQ=-.78125D-2*RQ*(VV-(4*K-1)**2)*(VV-(4*K+1)**2)/ 2296 & (K*(2*K+1)*X2) 2297 QX=QX+RQ 2298 ENDDO 2299 QX=.125D0*(VV-1D0)*QX/X 2300 XK=X-(.5D0*(J+V0)+.25D0)*PI 2301 A0=SQRT(RP2/X) 2302 CK=COS(XK) 2303 SK=SIN(XK) 2304 IF (J.EQ.0) THEN 2305 BJV0=A0*(PX*CK-QX*SK) 2306 ELSEIF (J.EQ.1) THEN 2307 BJV1=A0*(PX*CK-QX*SK) 2308 ENDIF 2309 ENDDO 2310 ENDIF 2311 BJ(0)=BJV0 2312 BJ(1)=BJV1 2313 IF (N.GE.2 .AND. N.LE.INT(.9*X)) THEN 2314 F0=BJV0 2315 F1=BJV1 2316 DO K=2,N 2317 F=2D0*(K+V0-1D0)/X*F1-F0 2318 BJ(K)=F 2319 F0=F1 2320 F1=F 2321 ENDDO 2322 ELSEIF (N.GE.2) THEN 2323 M=MSTA1(X,200) 2324 IF (M.LT.N) THEN 2325 N=M 2326 ELSE 2327 M=MSTA2(X,N,15) 2328 ENDIF 2329 F2=0D0 2330 F1=1D-100 2331 DO K=M,0,-1 2332 F=2D0*(V0+K+1D0)/X*F1-F2 2333 IF (K.LE.N) BJ(K)=F 2334 F2=F1 2335 F1=F 2336 ENDDO 2337 IF (ABS(BJV0).GT.ABS(BJV1)) THEN 2338 CS=BJV0/F 2339 ELSE 2340 CS=BJV1/F2 2341 ENDIF 2342 DO K=0,N 2343 BJ(K)=CS*BJ(K) 2344 ENDDO 2345 ENDIF 2346 VM=N+V0 2347 ENDIF 2348 END 2349C 2350 SUBROUTINE GAMMA(X,GA) 2351C 2352C ================================================== 2353C Purpose: Compute gamma function GA(x) 2354C Input : x --- Argument of GA(x) 2355C ( x is not equal to 0,-1,-2,... ) 2356C Output: GA --- GA(x) 2357C ================================================== 2358C 2359 IMPLICIT NONE 2360 INTEGER K,M,M1 2361 REAL*8 GA,GR,PI,R,X,Z 2362 PARAMETER(PI=3.141592653589793D0) 2363 REAL*8 G(26) 2364 DATA G/1D0, .5772156649015329D0, 2365 &-.6558780715202538D0, -.420026350340952D-1, 2366 &.1665386113822915D0, -.421977345555443D-1, 2367 &-.96219715278770D-2, .72189432466630D-2, 2368 &-.11651675918591D-2, -.2152416741149D-3, 2369 &.1280502823882D-3, -.201348547807D-4, 2370 &-.12504934821D-5, .11330272320D-5, 2371 &-.2056338417D-6, .61160950D-8, 2372 &.50020075D-8, -.11812746D-8, 2373 &.1043427D-9, .77823D-11, 2374 &-.36968D-11, .51D-12, 2375 &-.206D-13, -.54D-14, .14D-14, .1D-15/ 2376C 2377 IF (X.EQ.INT(X)) THEN 2378 IF (X.GT.0D0) THEN 2379 GA=1D0 2380 M1=X-1 2381 DO K=2,M1 2382 GA=GA*K 2383 ENDDO 2384 ELSE 2385 GA=1D300 2386 ENDIF 2387 ELSE 2388 IF (ABS(X).GT.1D0) THEN 2389 Z=ABS(X) 2390 M=INT(Z) 2391 R=1D0 2392 DO K=1,M 2393 R=R*(Z-K) 2394 ENDDO 2395 Z=Z-M 2396 ELSE 2397 Z=X 2398 ENDIF 2399 GR=G(26) 2400 DO K=25,1,-1 2401 GR=GR*Z+G(K) 2402 ENDDO 2403 GA=1D0/(GR*Z) 2404 IF (ABS(X).GT.1D0) THEN 2405 GA=GA*R 2406 IF (X.LT.0D0) GA=-PI/(X*GA*DSIN(PI*X)) 2407 ENDIF 2408 ENDIF 2409 END 2410C 2411 INTEGER FUNCTION MSTA1(X,MP) 2412C 2413C =================================================== 2414C Purpose: Determine the starting point for backward 2415C recurrence such that the magnitude of 2416C Jn(x) at that point is about 10^(-MP) 2417C Input : x --- Argument of Jn(x) 2418C MP --- Value of magnitude 2419C Output: MSTA1 --- Starting point 2420C =================================================== 2421C 2422 IMPLICIT NONE 2423 INTEGER IT,MP,N,N0,N1,NN 2424 REAL*8 A0,F,F0,F1,X 2425 REAL*8 ENVJ 2426 ENVJ(N,X)=.5D0*LOG10(6.28D0*N)-N*LOG10(1.36D0*X/N) 2427C 2428 A0=ABS(X) 2429 N0=INT(1.1D0*A0)+1 2430 F0=ENVJ(N0,A0)-MP 2431 N1=N0+5 2432 F1=ENVJ(N1,A0)-MP 2433 DO IT=1,20 2434 NN=N1-(N1-N0)/(1D0-F0/F1) 2435 F=ENVJ(NN,A0)-MP 2436 IF(ABS(NN-N1).LT.1) GOTO 20 2437 N0=N1 2438 F0=F1 2439 N1=NN 2440 F1=F 2441 ENDDO 244220 MSTA1=NN 2443 END 2444C 2445 INTEGER FUNCTION MSTA2(X,N,MP) 2446C 2447C =================================================== 2448C Purpose: Determine the starting point for backward 2449C recurrence such that all Jn(x) has MP 2450C significant digits 2451C Input : x --- Argument of Jn(x) 2452C n --- Order of Jn(x) 2453C MP --- Significant digit 2454C Output: MSTA2 --- Starting point 2455C =================================================== 2456C 2457 IMPLICIT NONE 2458 INTEGER IT,MP,N,N0,N1,NN 2459 REAL*8 A0,F,F0,F1,EJN,HMP,OBJ,X 2460 REAL*8 ENVJ 2461 ENVJ(N,X)=.5D0*LOG10(6.28D0*N)-N*LOG10(1.36D0*X/N) 2462C 2463 A0=ABS(X) 2464 HMP=.5D0*MP 2465 EJN=ENVJ(N,A0) 2466 IF (EJN.LE.HMP) THEN 2467 OBJ=MP 2468 N0=INT(1.1D0*A0) 2469 ELSE 2470 OBJ=HMP+EJN 2471 N0=N 2472 ENDIF 2473 F0=ENVJ(N0,A0)-OBJ 2474 N1=N0+5 2475 F1=ENVJ(N1,A0)-OBJ 2476 DO IT=1,20 2477 NN=N1-(N1-N0)/(1D0-F0/F1) 2478 F=ENVJ(NN,A0)-OBJ 2479 IF (ABS(NN-N1).LT.1) GOTO 20 2480 N0=N1 2481 F0=F1 2482 N1=NN 2483 F1=F 2484 ENDDO 248520 MSTA2=NN+10 2486 END 2487C 2488