NOTE

        THIS LIBRARY OF PROGRAMS FOR THE RECURSION METHOD WAS WRITTEN BY TOM

GODIN AND IS MADE AVAILABLE AT NO COST.  TOM AND I ONLY ASK THAT USERS

ACKNOWLEDGE THIS LIBRARY WHEN IT HAS MADE A SIGNIFICANT CONTRIBUTION

TO THEIR WORK, AND THAT THEY NOTIFY ME OF ANY ERRORS THEY DISCOVER.  I

STRONGLY ENCOURAGE USERS TO CONSULT WITH ME BEFORE EMBARKING ON LARGE PROJECTS

TO AVOID THE PITFALLS OTHERS HAVE ENCOUNTERED AND TO AVOID DUPLICATION OF

EFFORT.

        THESE PROGRAMS WERE ALSO PUBLISHED AS 'THE BLOCK RECURSION LIBRARY:

ACCURATE CALCULATIONS OF RESOLVENT SUBMATRICES USING THE BLOCK RECURSION

METHOD,' IN COMPUTER PHYSICS COMMUNICATIONS, VOLUME 64, PAGES 123-30, IN 1991.

        PLEASE ENJOY THE OPPORTUNITIES THESE PROGRAMS OFFER.





ROGER HAYDOCK

EUGENE, OREGON

2000 APRIL 6





0

0

0

0               OREGON QUANTUM CIRCUIT LIBRARY

+               ______________________________

0

0

0                               by

0                               T. Godin

                                Institute of Theoretical Science

                                University of Oregon

                                Eugene, OR 97403

                                USA

1        I.   Introduction

0             The Oregon Quantum Circuit Library (OQCL) is a series of

         programs and subroutines for numerically calculating certain

         electronic properties of devices in which quantum interference

         and/or tunneling effects are important. Each routine contains

         internal documentation explaining its function, and additional

         information is provided in this text.

0        II.  Physical Basis of the Calculation

0             The  routines in this library calculate  the  behavior of

         electrons (or other quantum mechanical particles) whose motion

         through a particular system is describable by a finite, dis-

         crete matrix form of an electronic Hamiltonian. This  Hamilto-

         nian  is  then given matrix elements with Hamiltonians repre-

         senting other regions of the circuit. By solving the correspon-

         ding eigenvalue problem,   we can obtain a quantum mechanical

         solution for the circuit as a whole; hence the name "Quantum

         Circuit Theory."

              A complete description of the theory behind such calcula-

         tions is beyond the scope of this introduction, but is pro-

         vided by the references. The calculations performed by the

         routines are described in Godin and Haydock, 1986. They are

         simple matrix generalizations of the recursion method, which

         is described in detail in Haydock, 1980.  The specific

         function of each routine is described in a separate section.

0        III.  Language Considerations

0             Routines in this library are written in a DEC-10 version

         of FORTRAN IV.  Most are written with the intent of avoiding

         extentions of the language that are peculiar to the DEC; how-

         ever, there are certain exceptions which users on other sys-

         tems will have to take into account. (These primarily relate

         to I/O.)

              In READ/WRITE statements, unit 5 refers to the terminal.

         Unit 6 refers to the lineprinter. Units 1-4 and 22-24 refer

         to disk files whose names are FORnn.DAT, where nn is the unit

         number refered to. Note that since these are default units on

         the DEC, no other disk commands (e.g., "OPEN") are used.  On

         other systems, such as IBM CMS, one would have to map the unit

         numbers to various devices.

              Other routines call library subroutines which are provided

         by the local computing facility. For example, a subroutine

         which generates a one-dimensional tight-binding lattice with

         random site energies calls an IMSL subroutine called GGUBS. If

         such a routine is not available another will have to substitu-

         ted. Also, some I/O routines which are DEC-specific have been

         deleted.

1             When a routine has one or more of these problems, it

         will be noted in the documentation.

              In updates and revisions to the original library, an

         effort has been made to use more standard I/O procedures and

         to segregate I/O to particular routines. It is hoped that

         this will facilitate conversion.

0        IV.  Disclaimer

0             There are certain limitations of this package that the

         user should be aware of.  Most of the routines were written

         for readability, not economy of code.  Also, optimization was

         not an overiding factor.

              Some sample programs have been included that utilize the

         routines in the library. These are not very sophisticated, are

         changed quite frequently and are not necessarily an integral

         part of the library. USE THEM AT YOUR OWN RISK. Such programs

         will be clearly marked as EXAMPLES.

              Hopefully, the OQCL subroutine library should be free of

         bugs. If there are any problems with them, errors, omissions

         or ambiguities in the documentation, compilation difficulties

         beyond normal dialect differences, I would appreciate hearing

         about it and will do what I can to help. If the user has any

         questions, comments or suggestions, please do not hesitate to

         contact me at the Institute of Theoretical Science, University

         of Oregon, Eugene OR 97403, telephone (503) 686 5222 (5204).

         BITNET ID: 55103 @ OREGON1.

             This library is distributed free of charge and may be used,

         copied or distributed at the discretion of the user, provided

         that appropriate citations are provided with the presentation

         of results obtained through its use.  The author disclaims any

         warrantee that the library be free of errors or for fitness

         for a particular purpose. No liability will be assumed for

         indirect, special, or consequntial damages.

1            LIBRARY FILES

0        All routines listed in this section which require user supplied

         routines or perform I/O will be preceded by an asterisk (*).

0        I.    BLOCKREC FORTRAN

0        BRECAL(N,M,NOVERM,A,B,THISA,THISB,THISP,NEXTB,D,ADJB,NEXTP,

                       LASTBP,NEXTBP)

         REAL A(N,N,NOVERM),B(N,N,NOVERM),THISA(N,N),THISB(N,N),

               THISP(M,N),NEXTB(M,N),D(M,N),ADJB(N,N),NEXTP(M,N),

                LASTBP(M,N),NEXTBP(M,N)

0        This subroutine takes the desired block dimension N, the

         Hamiltonian dimension M and their quotient NOVERM and tri-

         diagonalizes the Hamiltonian in NxN blocks. The  diagonal blocks

         are returned in the array A and the off diagonal blocks above

         the diagonal are returned in the array B. Each of B(i,j,1) is

         returned with a value of zero. All other arguments are dummies.

         (The off diagonal blocks below the diagonal are the adjoints

         of those above the diagonal since H is Hermitian.)

         BRECAL indirectly calls subroutine OFFDBL. The current version

         of this routine only allows tridiagonalization in 2x2 blocks.

         Modification of this routine will allow more general trans-

         formations.

0        FINDBL(LASTBP,NEXTBP,THISP,NEXTBP,THISP,NEXTP,THISA,

                     NEXTB,D,ADJB,N,M)

         REAL LASTBP(M,N),NEXTBP(M,N),THISP(M,N),NEXTP(M,N),

               THISA(N,N),NEXTB(N,N),D(M,N),ADJB(M,N)

0        This routine is the heart of the package. Its function is to

         actually perform one of the three-term recurrences, and it is

         called repeatedly by BRECAL until the recursion terminates,

         or is stopped. It takes the product of P(n-1) and the adjoint

         of B(n), and P(n) and computes A(n), P(n+1) and B(n+1). It

         calls in turn various specialized subroutines to perform these

         tasks. These values are input in, respectively, LASTBP, THISP

         and returned in THISA, THISP and NEXTB. The starting point

         for the next recursion, the product of the adjoint of B(n+1)

         and P(n), is returned in LASTBP.

1        *RDFHAM(M)

0        HAMOP(X,HX,N,M)

         REAL X(M,N),HX(M,N),H(20,20)

0        It is through these two routines that the user supplies the

         Hamiltonian to be transformed. RDFHAM reads a matrix of

         dimension M and passes it to HAMOP through a common block

         named HAM. The user may define RDFHAM to suit his system;

         a RDFHAM that works on a DEC 1091 and reads from a disk file

         named FULHAM.DAT is provided with the package.

         HAMOP takes a matrix X and operates with the Hamiltonian,

         returning the product in HX. This routine is called by HOP.

0

         HOP(X,Y,HX,A,N,M)

         REAL X(M,N),Y(M,N),HX(M,N),A(N,N)

0        This subroutine is a generalization of a similarly named

         routine in the Cambridge Recursion Library (Nex). It takes

         a matrix X and returns Y = H|X> - Y  and A = <X|H|X>.

0        DIAGBL(X,HX,A,N,M)

         REAL X(M,N),HX(M,N),A(N,N)

0        This routine, called by HOP, takes X and HX = H|X> and

         calculates A = <X|H|X>.

0        *OFFDBL(BP,N,M,P,B)

         REAL BP(M,N),P(M,N),B(N,N)

0        This subroutine performs the last step in the three-term

         recurrence. It takes the product of the next P, or next

         basis vector to be calculated, and the next off diagonal block

         B and separates them, returning them in variables of the same

         name.  The version of this routine supplied with the current

         package only works for N = 2. If another value is passed through

         this argument, an error message will be written on unit 5 and

         the program will halt execution. Generalization of this

         routine should be the only modification required to tridiagon-

         alize in arbitrary sized blocks.

1        ADJONT(B,ADJB,N)

         REAL B(N,N),ADJB(N,N)

0        The adjoint of matrix B is returned in ADJB.

0        VECSUB(X,Y,N,M)

         REAL X(M,N),Y(M,N)

0        The matrix X - Y is returned in Y.

0        BLMULT(A,X,AX,N,M)

         REAL A(N,N),X(M,N),AX(M,N)

0        This routine takes a block of the transformad Hamiltonian (A)

         and multiplies it by a vector represented in the basis of the

         untransformed Hamiltonian. The result, also expressed in the

         basis of the untransformed Hamiltonian, is returned in AX.

         This routine is useful, for example, in calculating products

         such as A(n)P(n) to subtract from HP(n).

0        INSERT(PIECE,WHOLE,N,L,I)

         REAL PIECE(N,N),WHOLE(N,N,L)

0        This routine takes a calculated matrix PIECE and loads it

         into an array of such matrices, so that WHOLE(j,k,I) =

         PIECE(j,k). This task is given to a routine to reduce clutter

         in BRECAL.

1        II.  TRALIB.FOR

0        This file contains I/O routines and routines for calculating

         transmission coefficients.

0        All of the follwing require the named common block

                  COMMON/TRANCO/A,B,N

         which contains the blocks of a 2x2 block tidiagonal Hamiltonian

         in A and B and the number of blocks in N. Ordinarily, the

         user need never see these locations or their contents; they

         should remain separate from the main program. The dimensions

         of A and B are (2,2,20) and (2,2,21).

0        The Routines:

0

0        *RDHAM

0        This subroutine reads the blocks of a 2x2 block-tridiagonal

         Hamiltonian from unit 1. Elements of the blocks must be

         unformatted, one block to a line, alternating between B's

         and A's. Data is terminated after allotted storage is filled

         or two lines of zero blocks (B=0,A=0) are encountered.

                  Sample:

                B1(1,1),(1,2),(2,1),(2,2)

                A1(1,1),  "  ,  "  ,  "

                B2...

                A2...

                0.,0.,0.,0.

                0.,0.,0.,0.

0        would enter 2 A's and 2 B's.

0

0        *CHART

0        This routine writes such a hamiltonian on unit 5 (DEC terminal)

0

0        *CHARTD

0        This routine has been removed from the package.

0

1        *COMPLEX FUNCTION TRANS(E)

         COMPLEX TRANS,E

0        If the hamiltonian read by RDHAM is a "conductor" which has

         two leads corresponding to the two diagonal sites in the first

         A block, then TRANS calculates the complex transmission

         coefficient of an electron incident on the conductor from one

         lead as a function of complex energy E.

         NOTE: TRANS writes error messages to unit 5.

0

0        COMPLEX FUNCTION TRANZ(Z)

         COMPLEX TRANZ,Z

0        TRANZ has been removed from the package.

0

0        III. Miscellaneous Routines

0

0        RANCHN

0        This subroutine generates a 1D tight-binding hamiltonian. At

         the user's instruction, either the site energies or the hopping

         or both will be random. the size of the fluctuations and the

         number of blocks are also specified. The user has the option of

         adding a constant potential to the system.

         NOTE: RANCHN reads instructions from unit 5 (terminal) and

         writes the result on unit 1 (disk) in a format that can be

         read by RDHAM.  It calls the IMSL subroutine GGUBS. If GGUBS is

         unavaillable a different routine will have to be substituted.

1        IV.   Example Programs

0        NOTE: These examples are just that. They are not very sophisti-

               cated and may contain bugs. They have lots of DEC I/O and

               one calls a plotting routine DRAW that was not included

               in the package.  PKSTAT usually gives reliable peak loca-

               tions but unreliable widths. USE THESE AT YOUR OWN RISK.

0        BLOCKREC FORTRAN

0        This program calls BRECAL and then writes the transformed

         Hamiltonian on unit 5. The untransformed Hamiltonian is read

         from a disk file called FULHAM.DAT by the subroutine RDFHAM

         which passes it via common to the subroutine HAMOP; the main

         program need never see it. The only function of the main

         program itself is to dimensionalize the dummy arguments passed

         to BRECAL and to output the transformed Hamiltonian.

0        MASTER FORTRAN

0        This program reads a block tridiagonal Hamiltonian from

         disk and calculates and plots the transmission coefficient

         of the system in the energy range the user specifies from

         unit 5 (terminal). It outputs its raw data on another disk

         file, which can be concanted with similar files and sent to

         a lineprinter. It also asks for two real serial numbers to  go
         with the data to distinguish the output from that of other

         runs. The user selects a regular, semilog, log or log-log plot.

0        PKSTAT(F,XDATA,N,PEAK,J)

         DIMENSION F(N),PEAK(50,4),XDATA(N),HALF

         REAL PEAK

0        This routine attempts to find peaks in data passed through

         F and XDATA. J gives the number of peaks found (up to 50).

         PEAK(I,1-4) gives height, location, width and width*height

         for each peak found. This is not a terribly sophisticated

         routine. Experience has shown that if the data has sufficient

         resolution, PKSTAT usually gives fairly reliable peak

         locations and heights; the values for the widths are not so

         good. The routine can also be "fooled" by unusual features

         in the data.

1        Recent Revisions to the Library

0           The current version of the library has many improvements,

         mainly in the block diagonalization routines.

            Serious errors in the previous block recursion program have

         been removed, and the whole package has been rewritten for

         clarity and ease of modification. With the generalization of

         one subroutine (OFFDBL) it will be cabable of performing the

         transformation in arbitrary sized blocks. These routines have

         been thoroughly tested and should perform satisfactorally.

            The transmission coefficient library has undergone minor

         modifications. The function TRANS has been modified to remove

         a minor error in the first version. The subroutines CHARTD

         and TRANZ have been removed on the assumption that the user

         can modify CHART or TRANS to fufill their functions as easily

         as he could adapt the supplied routines to the local FORTRAN

         I/O dialect.

            Adaptions to this package will probably continue as it is

         used and its strengths and weaknesses become apparent. Users

         are welcome to direct any comments along these lines to me.

0                                    Tom Godin

                                     November 10, 1986

0

0        REFERENCES

0        DEC System 10 FORTRAN Programmer's Reference Manual.

0        Godin, Haydock, "Quantum Circuit Theory," to be published in

                 J. Superlattices and Microstructures, 1986.

+                ____________________________________

0        Haydock, "The Recursive Solution of the Schroedinger Equation"

                 in Solid State Physics v.35, H. Ehrenreich, F.Seitz,

+                   ___________________

                 D. Turnbull eds., 216 (Academic Press, New York, 1980).

0        Nex, The Cambridge Recursion Library.

+             _______________________________













CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C

C    SUBROUTINE CHART, BY T. GODIN.

C

C     THIS SUBROUTINE PRINTS OUT THE BLOCKS OF A

C     LADDER-FORM HAMILTONIAN, I.E. THE NON-ZERO

C     BLOCKS OF A 2X2 BLOCK TRIDIAGONAL MATRIX.

C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C

      SUBROUTINE CHART

C

      INTEGER N

      REAL A(2,2,100), B(2,2,101)

      COMMON /TRANCO/ A,B,N

C

      WRITE(5,80)

80    FORMAT('0',11X,'A',27X,'B')

      DO 500 I=1, N

        WRITE(5,85) A(1,1,I),A(1,2,I),B(1,1,I+1),B(1,2,I+1)

        WRITE(5,85) A(2,1,I),A(2,2,I),B(2,1,I+1),B(2,2,I+1)

85      FORMAT(' ',2(2(E11.4,' '),4X))

        WRITE(5,86)

86      FORMAT('0')

500     CONTINUE

      RETURN

      END

C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C     RDHAM, a subroutine by T. Godin, U of O 11/25

C

C       This subroutine reads the blocks of a 2x2 block

C       tridiagonal Hamiltonian from a file called FOR01.DAT.

C

C

C    IMPORTANT---It is neccecarry to create a disk file

C      containing the input data which has the name FOR01.DAT

C      and contains this data in UNFORMATTED form. Such a

C      file is refered to in the program as unit 1.

C

C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

      SUBROUTINE RDHAM

C

       REAL A(2,2,100),B(2,2,101)

       INTEGER LIMIT

       COMMON /TRANCO/ A,B,N

       DATA LIMIT/100/

C

C  FIRST WE INPUT THE DATA, EACH LINE ALTERNATING A & B

C  SUBMATRICES.

C

       WRITE(5,11)

11     FORMAT(' Data being read from disk file FOR01.DAT....')

       N = 0

       B(1,1,1) = 1.

       B(1,2,1) = 0.

       B(2,1,1) = 0.

       B(2,2,1) = 1.

       A(1,1,1) = 0.

       A(1,2,1) = 0.

       A(2,1,1) = 0.

       A(2,2,1) = 0.

100    CONTINUE

         N = N+1

         READ(1,*)B(1,1,N+1),B(1,2,N+1),B(2,1,N+1),B(2,2,N+1)

         READ(1,*)A(1,1,N+1),A(1,2,N+1),A(2,1,N+1),A(2,2,N+1)

         DETER = B(1,1,N+1)*B(2,2,N+1) - B(1,2,N+1)*B(2,1,N+1)

         IF (DETER.NE.0.) GO TO 110

           WRITE(5,10)N,LIMIT

10         FORMAT(' ',I3,' blocks were read from disk. If the file',/,

     $      ' had more blocks that were not read,either your last',/,

     $      ' Bn had a zero eigenvalue or you exceeded storage',/,

     $      ' capacity. This routine stops when N = ',I3)

           GO TO 120

110        CONTINUE

         GO TO 100

120      CONTINUE

       RETURN

       END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C

C  TRANSMISSION COEFFICIENTS, A PROGRAM BY T. GODIN, U.OREGON, 7/85

C

C    THIS PROGRAM READS IN 2X2 BLOCKS OF A BLOCK-TRIDIAGONAL MATRIX

C     AND COMPUTES THE TRANSMISSION COEFFICIENTS OF ELECTRONS OF

C      GIVEN ENERGY IMPINGING ON A POTENTIAL REPRESENTED BY SUCH

C      A MATRIX

C

C    This function is the "guts" of the old TRANS program modified

C   to accept complex energies. It gives t[E] only, and lets the

C   main decide what E's {which might lie in complex plane}.

C

C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

      COMPLEX FUNCTION TRANS(ENERGY)

C

       REAL A(2,2,100),B(2,2,101),BINV(2,2),DETER,ADJB(2,2)

       COMPLEX X1(2,2),X2(2,2),Y1(2,2),X

       COMPLEX Y2(2,2),ELESSA(2,2),ENERGY,Z1(2,2),Z2(2,2)

       COMPLEX C(2),D(2,2),DINV(2,2),CDET,CTRANS

       INTEGER N

       COMMON /TRANCO/ A,B,N

C

C  WE MUST FIRST SET THE INITIAL CONDITIONS...

C

         DO 200 I=1,2

           DO 210 J=1,2

             X1(I,J)=(0.,0.)

             Y1(I,J)=(0.,0.)

             X2(I,J)=(0.,0.)

             Y2(I,J)=(0.,0.)

210          CONTINUE

           X2(I,I)=(1.,0.)

           Y1(I,I)=(1.,0.)

200        CONTINUE

C

C  NOW WE WILL BE USING THESE VALUES TO PERFORM THE RECCURENCE

C  RELATION FOR THE WAVEFUNCTIONS. THE RELATION TERMINATES WHEN

C  WE GET A B WITH A ZERO EIGENVALUE, AT WHICH POINT THE FINAL

C  RECCURENCE (DIFFERENT IN FORM) IS CALCULATED.

C

         DO 300 I=1,N

           DETER = B(1,1,I+1)*B(2,2,I+1) - B(1,2,I+1)*B(2,1,I+1)

           IF (DETER.EQ.0.) GO TO 315

           BINV(1,1) =  B(2,2,I+1)/DETER

           BINV(1,2) = -B(1,2,I+1)/DETER

           BINV(2,1) = -B(2,1,I+1)/DETER

           BINV(2,2) =  B(1,1,I+1)/DETER

315        CONTINUE

C

C  NOW WE SET UP THE MEMORY LOCATIONS NEEDED FOR THE CALCULATION.

C  ELESSA(I,J) CONTAINS THE MATRIX (EI-A), THE Z'S ARE USED AS TEMPORARY

C  STORAGE DURING THE CALCULATION, AND THEN THE FINAL VALUES

C  ARE STORED THERE AT THE END.

C

           DO 320 J=1,2

             DO 330 K=1,2

               ELESSA(J,K) = -A(J,K,I)*(1.,0.)

               Z1(J,K) = (0.,0.)

               Z2(J,K) = (0.,0.)

330            CONTINUE

             ELESSA(J,J) = ELESSA(J,J)+ENERGY

320          CONTINUE

C

           DO 333 J=1,2

             DO 334 K=1,2

               ADJB(J,K) = B(K,J,I)

334            CONTINUE

333          CONTINUE

C

           DO 340 J=1,2

             DO 350 K=1,2

               DO 355 L=1,2

                Z1(J,K)=Z1(J,K)+ELESSA(J,L)*X2(L,K)-ADJB(J,L)*X1(L,K)

                Z2(J,K)=Z2(J,K)+ELESSA(J,L)*Y2(L,K)-ADJB(J,L)*Y1(L,K)

355              CONTINUE

350            CONTINUE

340          CONTINUE

C

C  CHECK TO SEE IF THIS IS THE LAST RECCURSION; IF SO IT IS NOT

C  NECCESSARY TO MULTIPLY BY B INVERSE.

C

           IF (DETER.EQ.0.) GO TO 399

           DO 360 J=1,2

             DO 370 K=1,2

               X1(J,K)=X2(J,K)

               Y1(J,K)=Y2(J,K)

               X2(J,K)=(0.,0.)

               Y2(J,K)=(0.,0.)

               DO 380 L=1,2

                 X2(J,K)=X2(J,K)+BINV(J,L)*Z1(L,K)

                 Y2(J,K)=Y2(J,K)+BINV(J,L)*Z2(L,K)

380              CONTINUE

370            CONTINUE

360          CONTINUE

399        CONTINUE

300        CONTINUE

C

C  NOW WE PERFORM THE LAST BIT OF MATRIX ALGEBRA TO GET t ,

C    AND SEND IT BACK TO THE MAIN PROGRAM.

C

         X = CSQRT( (4.,0.) - ENERGY*ENERGY)

         DO 400 J=1,2

           C(J)=Z1(J,1)*(ENERGY + (0.,1.)*X)/2. + Z2(J,1)

           DO 410 K=1,2

             D(J,K)=Z1(J,K)*(ENERGY + (0.,-1.)*X)/2. + Z2(J,K)

410          CONTINUE

400        CONTINUE

         CDET = D(1,1)*D(2,2) - D(1,2)*D(2,1)

         IF(CDET.NE.(0.,0.)) GO TO 420

           WRITE(5,411) ENERGY

411        FORMAT(' ERROR: BOUNDARY CONDITIONS CANNOT BE',/,

     $          ' MATCHED AT E = ',F4.2)

           TRANS = (0.,0.)

           GO TO 430

420        CONTINUE

         DINV(2,1)=-D(2,1)/CDET

         DINV(2,2)=+D(1,1)/CDET

         TRANS = -DINV(2,1)*C(1) - DINV(2,2)*C(2)

430      CONTINUE

       RETURN

       END









      SUBROUTINE ADJONT (B,ADJB,N)

      REAL B(N,N),ADJB(N,N)

      DO 10 I=1, N

        DO 20 J=1,N

          ADJB(I,J) = B(J,I)

20        CONTINUE

10      CONTINUE

      RETURN

      END

C

      SUBROUTINE OFFDBL(PROD,N,M,P,B)

      REAL PROD(M,N),P(M,N),B(2,2)

C

      IF (N.EQ.2) GO TO 100

        WRITE(5,10)

10      FORMAT(' ***ERROR--Global module OFFDBL only works for 2x2

     + blocks.')

        STOP

100     CONTINUE

C

C   WE USE THE CONVENTION THAT B(1,2) = 0. THIS FIXES B.

C

      B(1,2) = 0.

C

      B(1,1) = 0.

      DO 200 I = 1, M

        B(1,1) = B(1,1) + PROD(I,1)*PROD(I,1)

200     CONTINUE

      B(1,1) = SQRT(B(1,1))

      IF (B(1,1).EQ.0.) GO TO 210

      DO 210 I=1, M

        P(I,1) = PROD(I,1)/B(1,1)

210     CONTINUE

C

      B(2,1) = 0.
      DO 300 I=1, M

        B(2,1) = B(2,1) + P(I,1)*PROD(I,2)

300     CONTINUE

C

      B(2,2) = 0.

      DO 400 I=1, M

        P(I,2) = PROD(I,2) - B(2,1)*P(I,1)

        B(2,2) = B(2,2) + P(I,2)*P(I,2)

400     CONTINUE

      B(2,2) = SQRT(B(2,2))

      IF (B(2,2).EQ.0.) GO TO 410

      DO 410 I=1, M

        P(I,2) = P(I,2)/B(2,2)

410     CONTINUE

C

      RETURN

      END

C

      SUBROUTINE VECSUB (X,Y,N,M)

      REAL X(M,N), Y(M,N)

      DO 10 I=1, M

        DO 11 J=1, N

          Y(I,J) = X(I,J) - Y(I,J)

11        CONTINUE

10      CONTINUE

      RETURN

      END

C

      SUBROUTINE DIAGBL(X,Y,A,N,M)

      REAL X(M,N),Y(M,N),A(N,N)

C

      DO 100 I=1,N

        DO 110 J=1,N

          A(I,J) = 0.

110       CONTINUE

100     CONTINUE

C

      DO 200 I=1,N

        DO 210 J=1, N

          DO 220 K=1, M

            A(I,J) = A(I,J) + X(K,I)*Y(K,J)

220         CONTINUE

210       CONTINUE

200     CONTINUE

      RETURN

      END

C

      SUBROUTINE HOP (X,Y,HX,A,N,M)

      REAL X(M,N), Y(M,N), A(N,N),HX(M,N)

C

      CALL HAMOP(X,HX,N,M)

C

      DO 100 I=1,M

        DO 110 J=1,N

          Y(I,J) = HX(I,J) - Y(I,J)

110       CONTINUE

100     CONTINUE

C

      CALL DIAGBL(X,HX,A,N,M)

C

      RETURN

      END

C

      SUBROUTINE HAMOP(X,HX,N,M)

      REAL X(M,N),HX(M,N),H(20,20)

      COMMON/HAM/ H

      K = 20
C

      DO 100 I=1,M

        DO 110 J=1,N

          HX(I,J) = 0.

110       CONTINUE

100     CONTINUE

C

      DO 200 I=1,M

        DO 210 J=1, N

          DO 220 L=1, K

            HX(I,J) = HX(I,J) + H(I,L)*X(L,J)

220         CONTINUE

210       CONTINUE

200     CONTINUE

C

      RETURN

      END

C

      SUBROUTINE RDFHAM(M)

      REAL H(20,20)

      COMMON /HAM/ H

      IF (M.LE.20) GO TO 27

        WRITE (5,28)

28      FORMAT(' ***Hamiltonian on disk too large for allotted space.')

        STOP

27      CONTINUE

C

      OPEN (UNIT=30,DEVICE='DSK',FILE='FULHAM.DAT')

      DO 10 I=1,M

        READ(30,*)(H(I,J),J=1,M)

10      CONTINUE

      CLOSE(UNIT=30)

      RETURN

      END

C

      SUBROUTINE BLMULT(A,X,AX,N,M)

      REAL A(N,N),X(M,N),AX(M,N)

C

      DO 20 I=1,M

        DO 21 J=1,N

          AX(I,J) = 0.

21        CONTINUE

20      CONTINUE

C

      DO 10 I=1, N

        DO 11 J=1, N

          DO 12 K=1, M

            AX(K,I) = AX(K,I) + A(I,J)*X(K,J)

12          CONTINUE

11        CONTINUE

10      CONTINUE

      RETURN

      END

C

      SUBROUTINE INSERT (PIECE,WHOLE,N,L,I)

      REAL PIECE(N,N),WHOLE(N,N,L)

C

      DO 10 J=1, N

        DO 11 K=1, N

          WHOLE (J,K,I) = PIECE(J,K)

11        CONTINUE

10      CONTINUE

      RETURN

      END

C

      SUBROUTINE FINDBL(LASTBP,NEXTBP,THISP,NEXTP, THISA,NEXTB,D,ADJB,

     +                          N,M)

      REAL LASTBP(M,N),NEXTBP(M,N),THISP(M,N),NEXTP(M,N),THISA(N,N),

     +        NEXTB(N,N),D(M,N),ADJB(N,N)

C

C*******1. Get HPn - BnPn-1 into LASTBP, find THISA.*****************

      CALL HOP(THISP,LASTBP,D,THISA,N,M)

C*******2. Find AnPn.             ********************************

      CALL BLMULT(THISA,THISP,NEXTBP,N,M)

C*******3. Find NEXTBP = HPn - BnPn-1 - AnPn  =  Bn+1Pn+1. *******

      CALL VECSUB(LASTBP,NEXTBP,N,M)

C*******4. Separate the product Bn+1Pn+1 into Bn+1,  Pn+1. ********

      CALL OFFDBL(NEXTBP,N,M,NEXTP,NEXTB)

C*******5. Finally, set up Bn+1Pn  and Pn in appropriate variables

C                       for the next recursion.*************************

      CALL ADJONT(NEXTB,ADJB,N)

      CALL BLMULT(ADJB,THISP,LASTBP,N,M)

        write(5,102)((thisp(ii,jj),jj=1,n),ii=1,m)

102     format(' p:',8(/,' ',2(e13.2)))

      DO 100 I=1,M

        DO 110 J=1,N

          THISP(I,J) = NEXTP(I,J)

110       CONTINUE

100     CONTINUE

      RETURN

      END

C

      SUBROUTINE BRECAL(N,M,N OVER M, A,B,THISA,THISB,THISP,NEXTB,

     +                 D,ADJB,NEXTP,LASTBP,NEXTBP)

      REAL A(N,N,N OVER M),B(N,N,N OVER M),THISA(N,N),THISB(N,N),

     +    THISP(M,N),NEXTP(M,N),NEXTB(N,N),LASTBP(M,N),NEXTBP(M,N),

     +                                          D(M,N),ADJB(N,N)

C

      CALL RDFHAM(M)

C

C   FIRST WE SET UP INITIAL CONDITIONS AND PERFORM THE FIRST RECURSION.

C

      I = 1

      DO 100 J=1, N

        DO 110 K=1, M

          LASTBP(K,J) = 0.

          THISP(K,J) = 0.

110       CONTINUE

        THISP(J,J) = 1.

100     CONTINUE

C

      CALL FINDBL(LASTBP,NEXTBP,THISP,NEXTP,THISA,NEXTB,D,ADJB,N,M)

      CALL INSERT(THISA,A,N,NOVERM,I)

      CALL INSERT(NEXTB,B,N,NOVERM,I+1)

C

C   NOW THE MAIN LOOP WILL CALCULATE SUCCESSIVE LEVELS. THE INITIAL

C     CONDITIONS FOR EACH ITERATIONS ARE TAKEN CARE OF BY THE SUBROUTINES.

C

      DO 200 I=2, N OVER M - 1

        J = I

        CALL FINDBL(LASTBP,NEXTBP,THISP,NEXTP,THISA,NEXTB,D,ADJB,N,M)

        CALL INSERT(THISA,A,N,NOVERM,J)

        CALL INSERT(NEXTB,B,N,NOVERM,J+1)

200     CONTINUE

C

C   NOW THE LAST RECURSION. WE DON'T NEED TO STORE B HERE.

C

      I = NOVERM

      CALL FINDBL(LASTBP,NEXTBP,THISP,NEXTP,THISA,NEXTB,D,ADJB,N,M)

      CALL INSERT(THISA,A,N,NOVERM,I)

      RETURN

      END

ÿ