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 = .
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 = .
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
ÿ