MultiPhysicsVault/.raw/FiniteElementProcedures/FiniteElementProcedures_099.md

<!-- source-page: 981 -->

![](images/page-981_763fe56350c9f90ac455b11e1132c0745e2a2fa873c09fcf8bab94adddf4e3f7.jpg)

<details>
<summary>text_image</summary>

Bound on λ₁
Bound on λ₄ and λ₅
μ
λ
λ₂⁽ʲ⁺¹⁾
Radius 0.01 λ₄⁽ʲ⁺¹⁾
</details>

Figure 11.5 Bounds on eigenvalues to apply Sturm sequence check, p = 6

erty of the characteristic polynomials of problems $\mathbf{K}\boldsymbol{\phi} = \lambda \mathbf{M}\boldsymbol{\phi}$ and $\mathbf{K}^{(r)}\boldsymbol{\phi}^{(r)} = \lambda^{(r)}\mathbf{M}^{(r)}\boldsymbol{\phi}^{(r)}$ at a shift $\mu$ , where $\mu$ is just to the right of the calculated value for $\lambda_p$ (see Fig. 11.5). The Sturm sequence property yields that in the Gauss factorization of $\mathbf{K} - \mu \mathbf{M}$ into $\mathbf{LDL}^T$ , the number of negative elements in $\mathbf{D}$ is equal to the number of eigenvalues smaller than $\mu$ . Hence, in the case considered, we should have $p$ negative elements in $\mathbf{D}$ . However, in order to apply the Sturm sequence check, a meaningful $\mu$ must be used that takes account of the fact that we have obtained only approximations to the exact eigenvalues of the problem $\mathbf{K}\boldsymbol{\phi} = \lambda \mathbf{M}\boldsymbol{\phi}$ . Let $l$ be the last iteration, so that the calculated eigenvalues are $\lambda_1^{(l+1)}$ , $\lambda_2^{(l+1)}$ , $\ldots$ , $\lambda_p^{(l+1)}$ . Since (11.156) is satisfied, we can use

$$
0. 9 9 \lambda_ {i} ^ {(l + 1)} \leq \lambda_ {i} <   1. 0 1 \lambda_ {i} ^ {(l + 1)} \tag {11.157}
$$

or tighter bounds based on the actual accuracy reached in (11.156). The relation in (11.157) can then be used to establish bounds on all exact eigenvalues, and hence a realistic Sturm sequence check can be applied.

# 11.6.5 Implementation of the Subspace Iteration Method

The equations of subspace iteration have been presented in (11.151) to (11.155). However, in actual implementation, the solution can be performed more effectively as summarized in Table 11.3, which for p small also gives the corresponding number of operations used.

The solution method is presented in a compact manner in the computer program SSPACE. This program provides only an implementation of the basic steps of the subspace iteration method described above without including acceleration techniques that are important in practice. One important aspect of the method is its relative simplicity when compared with other solution techniques, and this simplicity is also reflected in the program SSPACE.

Subroutine SSPACE. Program SSPACE is an implementation of the basic subspace iteration method presented above for the solution of the smallest eigenvalues and corresponding eigenvectors of the generalized eigenproblem $K\phi = \lambda M\phi$ . The argument variables and use of the subroutine are defined using comment lines in the program.

<!-- source-page: 982 -->

TABLE 11.3 Summary of subspace iteration solution, $p \leq  {20}$ 

<table><tr><td rowspan="2">Operation</td><td rowspan="2">Calculation</td><td colspan="2">Number of operations</td><td rowspan="2">Required storage</td></tr><tr><td> $m = m_{\mathrm{K}} = m_{\mathrm{M}}$ </td><td> $m = m_{\mathrm{K}}, m_{\mathrm{M}} = 0$ </td></tr><tr><td>Factorization of K</td><td> $\mathbf{K} = \mathbf{LDL}^{T}$ </td><td> $\frac{1}{2}nm^{2} + \frac{3}{2}nm$ </td><td> $\frac{1}{2}nm^{2} + \frac{3}{2}nm$ </td><td></td></tr><tr><td rowspan="6">Subspace iteration</td><td> $\mathbf{K} \overline{\mathbf{X}}_{k+1} = \mathbf{Y}_{k}$ </td><td> $nq(2m + 1)$ </td><td> $nq(2m + 1)$ </td><td></td></tr><tr><td> $\mathbf{K}_{k+1} = \overline{\mathbf{X}}_{l+1}^{T} \mathbf{Y}_{k}$ </td><td> $\frac{1}{2}nq(q + 1)$ </td><td> $\frac{1}{2}nq(q + 1)$ </td><td></td></tr><tr><td> $\overline{\mathbf{Y}}_{k+1} = \mathbf{M} \overline{\mathbf{X}}_{k+1}$ </td><td> $nq(2m + 1)$ </td><td> $nq$ </td><td></td></tr><tr><td> $\mathbf{M}_{k+1} = \overline{\mathbf{X}}_{l+1}^{T} \overline{\mathbf{Y}}_{k+1}$ </td><td> $\frac{1}{2}nq(q + 1)$ </td><td> $\frac{1}{2}nq(q + 1)$ </td><td rowspan="3">Algorithm is effectively implemented as out-of-core solver</td></tr><tr><td> $\mathbf{K}_{k+1} \mathbf{Q}_{k+1} = \mathbf{M}_{k+1} \mathbf{Q}_{k+1} \mathbf{A}_{k+1}$ </td><td></td><td> $o(q^{3})$  neglected</td></tr><tr><td> $\mathbf{Y}_{k+1} = \overline{\mathbf{Y}}_{k+1} \mathbf{Q}_{k+1}$ </td><td> $nq^{2}$ </td><td> $nq^{2}$ </td></tr><tr><td rowspan="2">Sturm sequence check</td><td> $\overline{\mathbf{K}} = \mathbf{K} - \mu \mathbf{M}$ </td><td> $n(m + 1)$ </td><td> $n$ </td><td></td></tr><tr><td> $\overline{\mathbf{K}} = \mathbf{LDL}^{T}$ </td><td> $\frac{1}{2}nm^{2} + \frac{3}{2}nm$ </td><td> $\frac{1}{2}nm^{2} + \frac{3}{2}nm$ </td><td></td></tr><tr><td>Error* measures</td><td> $\frac{\|\mathbf{K}\phi_{i}^{(l+1)} - \lambda_{i}^{(l+1)} \mathbf{M}\phi_{i}^{(l+1)}\|_{2}}{\|\mathbf{K}\phi_{i}^{(l+1)}\|_{2}}$ </td><td> $4nm + 5n$ </td><td> $2nm + 5n$ </td><td></td></tr><tr><td colspan="2">Total for solution of  $p$  lowest eigenvalues and associated eigenvectors, assuming that ten iterations are required and  $q = \min \{2p, p + 8\}$ </td><td> $nm^{2} + nm(4 + 4p) + 5np + 20nq(2m + q + \frac{3}{2})$ </td><td> $nm^{2} + nm(3 + 2p) + 5np + 20nq(m + q + \frac{3}{2})$ </td><td></td></tr></table>

\*The error measures are not needed but may be of interest.

<!-- source-page: 983 -->

SUBROUTINE SSPACE (A, B, MAXA, R, EIGV, TT, W, AR, BR, VEC, D, RTOLV, BUP, BLO, SSP000011 BUPC, NN, NNM, NWK, NWM, NROOT, RTOL, NC, NNC, NITEM, IFSS, IFPR, NSTIF, IOUT) SSP00002

# PROGRAM

TO SOLVE FOR THE SMALLEST EIGENVALUES-- ASSUMED .GT. 0 -- AND CORRESPONDING EIGENVECTORS IN THE GENERALIZED EIGENPROBLEM USING THE SUBSPACE ITERATION METHOD

# -- INPUT VARIABLES --

<table><tr><td>A(NWK)</td><td>= STIFFNESS MATRIX IN COMPACTED FORM (ASSUMED POSITIVE DEFINITE)</td></tr><tr><td>B(NWM)</td><td>= MASS MATRIX IN COMPACTED FORM</td></tr><tr><td>MAXA(NNM)</td><td>= VECTOR CONTAINING ADDRESSES OF DIAGONAL ELEMENTS OF STIFFNESS MATRIX A</td></tr><tr><td>R(NN,NC)</td><td>= STORAGE FOR EIGENVECTORS</td></tr><tr><td>EIGV(NC)</td><td>= STORAGE FOR EIGENVALUES</td></tr><tr><td>TT(NN)</td><td>= WORKING VECTOR</td></tr><tr><td>W(NN)</td><td>= WORKING VECTOR</td></tr><tr><td>AR(NNC)</td><td>= WORKING MATRIX STORING PROJECTION OF K</td></tr><tr><td>BR(NNC)</td><td>= WORKING MATRIX STORING PROJECTION OF M</td></tr><tr><td>VEC(NC,NC)</td><td>= WORKING MATRIX</td></tr><tr><td>D(NC)</td><td>= WORKING VECTOR</td></tr><tr><td>RTOLV(NC)</td><td>= WORKING VECTOR</td></tr><tr><td>BUP(NC)</td><td>= WORKING VECTOR</td></tr><tr><td>BLO(NC)</td><td>= WORKING VECTOR</td></tr><tr><td>BUPC(NC)</td><td>= WORKING VECTOR</td></tr><tr><td>NN</td><td>= ORDER OF STIFFNESS AND MASS MATRICES</td></tr><tr><td>NNM</td><td>= NN + 1</td></tr><tr><td>NWK</td><td>= NUMBER OF ELEMENTS BELOW SKYLINE OF STIFFNESS MATRIX</td></tr><tr><td>NWM</td><td>= NUMBER OF ELEMENTS BELOW SKYLINE OF MASS MATRIX</td></tr></table>

I. E. NWM=NWK FOR CONSISTENT MASS MATRIX
NWM=NN FOR LUMPED MASS MATRIX

<table><tr><td>NROOT</td><td>= NUMBER OF REQUIRED EIGENVALUES AND EIGENVECTORS</td></tr><tr><td>RTOL</td><td>= CONVERGENCE TOLERANCE ON EIGENVALUES( 1.E-06 OR SMALLER )</td></tr></table>

<table><tr><td>NC</td><td>= NUMBER OF ITERATION VECTORS USED(USUALLY SET TO MIN(2*NROOT, NROOT+8), BUT NCCANNOT BE LARGER THAN THE NUMBER OF MASSDEGREES OF FREEDOM)</td></tr></table>

<table><tr><td>NNC</td><td>= NC*(NC+1)/2 DIMENSION OF STORAGE VECTORS AR,BR</td></tr></table>

<table><tr><td>NITEM</td><td>= MAXIMUM NUMBER OF SUBSPACE ITERATIONS PERMITTED(USUALLY SET TO 16)THE PARAMETERS NC AND/OR NITEM MUST BEINCREASED IF A SOLUTION HAS NOT CONVERGED</td></tr><tr><td>IFSS</td><td>= FLAG FOR STURM SEQUENCE CHECKEQ.0 NO CHECKEQ.1 CHECK</td></tr></table>

<table><tr><td>IFPR</td><td>= FLAG FOR PRINTING DURING ITERATION</td></tr><tr><td></td><td>EQ.0 NO PRINTING</td></tr><tr><td></td><td>EQ.1 PRINT</td></tr></table>

<table><tr><td>NSTIF</td><td>= SCRATCH FILE</td></tr><tr><td>IOUT</td><td>= UNIT USED FOR OUTPUT</td></tr></table>

# -- OUTPUT --

EIGV(NROOT) = EIGENVALUES
R(NN, NROOT) = EIGENVECTORS

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON CRAY EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM MACHINES, ENGINEERING WORKSTATIONS AND PCS. DEACTIVATE ABOVE LINE FOR SINGLE PRECISION ARITHMETIC.

INTEGER MAXA(NNM)

DIMENSION A(NWK), B(NWM), R(NN, NC), TT(NN), W(NN), EIGV(NC),

<!-- source-page: 984 -->

```csv
1 D(NC), VEC(NC, NC), AR(NNC), BR(NNC), RTOLV(NC), BUP(NC), SSP00071
2 BLO(NC), BUPC(NC) SSP00072
C
C SET TOLERANCE FOR JACOBI ITERATION SSP00073
C TOLJ=1.0D-12 SSP00074
C
C INITIALIZATION SSP00075
C
ICONV=0 SSP00076
NSCH=0 SSP00077
NSMAX=12 SSP00078
N1=NC + 1 SSP00079
NC1=NC - 1 SSP00080
REWIND NSTIF SSP00081
WRITE (NSTIF) A SSP00082
DO 2 I=1, NC SSP00083
DO 2 I=1, NC SSP00084
2 D(I)=0. SSP00085
C
C ESTABLISH STARTING ITERATION VECTORS SSP00086
C
ND=NN/NC SSP00087
IF (NWM.GT.NN) GO TO 4 SSP00088
J=0 SSP00089
DO 6 I=1, NN SSP00090
II=MAXA(I) SSP00091
R(I,1)=B(I) SSP00092
IF (B(I).GT.0) J=J + 1 SSP00093
6 W(I)=B(I)/A(II) SSP00094
IF (NC.LE.J) GO TO 16 SSP00095
WRITE (IOUT,1007) SSP00096
GO TO 800 SSP00097
4 DO 10 I=1, NN SSP00098
II=MAXA(I) SSP00099
R(I,1)=B(II) SSP00100
10 W(I)=B(II)/A(II) SSP00101
16 DO 20 J=2, NC SSP00102
DO 20 I=1, NN SSP00103
20 R(I,J)=0. SSP00104
C
L=NN - ND SSP00105
DO 30 J=2, NC SSP00106
RT=0. SSP00107
DO 40 I=1, L SSP00108
IF (W(I).LT.RT) GO TO 40 SSP00109
RT=W(I) SSP00110
IJ=I SSP00111
40 CONTINUE SSP00112
DO 50 I=L, NN SSP00113
IF (W(I).LE.RT) GO TO 50 SSP00114
RT=W(I) SSP00115
IJ=I SSP00116
50 CONTINUE SSP00117
TT(J)=FLOAT(IJ) SSP00118
W(IJ)=0. SSP00119
L=L - ND SSP00120
30 R(IJ,J)=1. SSP00121
C
WRITE (IOUT,1008) SSP00122
WRITE (IOUT,1002) (TT(J),J=2, NC) SSP00123
A RANDOM VECTOR IS ADDED TO THE LAST VECTOR SSP00124
C
PI=3.141592654D0 SSP00125
XX=0.5D0 SSP00126
DO 60 K=1, NN SSP00127
XX=(PI + XX)**5 SSP00128
IX=INT(XX) SSP00129
XX=XX - FLOAT(IX) SSP00130
60 R(K,NC)=R(K,NC) + XX SSP00131
C 
```

<!-- source-page: 985 -->

```csv
C FACTORIZE MATRIX A INTO (L)*(D)*(L(T))
C
ISH=0
CALL DECOMP (A,MAXA,NN,ISH,IOUT)
C
C --- - START OF ITERATION LOOP
C
NITE=0
TOLJ2=1.0D-24
100 NITE=NITE + 1
IF (IFPR.EQ.0) GO TO 90
WRITE (IOUT,1010) NITE
C
CALCULATE THE PROJECTIONS OF A AND B
C
90 IJ=0
DO 110 J=1,NC
DO 120 K=1,NN
120 TT(K)=R(K,J)
CALL REDBAK (A,TT,MAXA,NN)
DO 130 I=J,NC
ART=0.
DO 140 K=1,NN
140 ART=ART + R(K,I)*TT(K)
IJ=IJ + 1
130 AR(IJ)=ART
DO 150 K=1,NN
150 R(K,J)=TT(K)
110 CONTINUE
IJ=0
DO 160 J=1,NC
CALL MULT (TT,B,R(1,J),MAXA,NN,NWM)
DO 180 I=J,NC
BRT=0.
DO 190 K=1,NN
190 BRT=BRT + R(K,I)*TT(K)
IJ=IJ + 1
180 BR(IJ)=BRT
IF (ICONV.GT.0) GO TO 160
DO 200 K=1,NN
200 R(K,J)=TT(K)
160 CONTINUE
C
SOLVE FOR EIGENSYSTEM OF SUBSPACE OPERATORS
C
IF (IFPR.EQ.0) GO TO 320
IND=1
210 WRITE (IOUT,1020)
II=1
DO 300 I=1,NC
ITEMP=II + NC - I
WRITE (IOUT,1005) (AR(J),J=II,ITEMP)
300 II=II + N1 - I
WRITE (IOUT,1030)
II=1
DO 310 I=1,NC
ITEMP=II + NC - I
WRITE (IOUT,1005) (BR(J),J=II,ITEMP)
310 II=II + N1 - I
IF (IND.EQ.2) GO TO 350
C
320 CALL JACOBI (AR,BR,VEC,EIGV,W,NC,NNC,TOLJ,NSMAX,IFPR,IOUT)
C
IF (IFPR.EQ.0) GO TO 350
WRITE (IOUT,1040)
IND=2
GO TO 210
C
ARRANGE EIGENVALUES IN ASCENDING ORDER
C
SSP00141
SSP00142
SSP00143
SSP00144
SSP00145
SSP00146
SSP00147
SSP00148
SSP00149
SSP00150
SSP00151
SSP00152
SSP00153
SSP00154
SSP00155
SSP00156
SSP00157
SSP00158
SSP00159
SSP00160
SSP00161
SSP00162
SSP00163
SSP00164
SSP00165
SSP00166
SSP00167
SSP00168
SSP00169
SSP00170
SSP00171
SSP00172
SSP00173
SSP00174
SSP00175
SSP00176
SSP00177
SSP00178
SSP00179
SSP00180
SSP00181
SSP00182
SSP00183
SSP00184
SSP00185
SSP00186
SSP00187
SSP00188
SSP00189
SSP00190
SSP00191
SSP00192
SSP00193
SSP00194
SSP00195
SSP00196
SSP00197
SSP00198
SSP00199
SSP00200
SSP00201
SSP00202
SSP00203
SSP00204
SSP00205
SSP00206
SSP00207
SSP00208
SSP00209
SSP00210 
```

<!-- source-page: 986 -->

```csv
350 IS=0
    II=1
    DO 360 I=1,NC1
    ITEMP=II + N1 - I
    IF (EIGV(I+1).GE.EIGV(I)) GO TO 360
    IS=IS + 1
    EIGVT=EIGV(I+1)
    EIGV(I+1)=EIGV(I)
    EIGV(I)=EIGVT
    BT=BR(ITEMP)
    BR(ITEMP)=BR(II)
    BR(II)=BT
    DO 370 K=1,NC
    RT=VEC(K,I+1)
    VEC(K,I+1)=VEC(K,I)
370 VEC(K,I)=RT
360 II=ITEMP
    IF (IS.GT.0) GO TO 350
    IF (IFPR.EQ.0) GO TO 375
    WRITE (IOUT,1035)
    WRITE (IOUT,1006) (EIGV(I),I=1,NC)
C
    CALCULATE B TIMES APPROXIMATE EIGENVECTORS (ICONV.EQ.0)
C    OR    FINAL EIGENVECTOR APPROXIMATIONS (ICONV.GT.0)
C
375 DO 420 I=1,NN
    DO 422 J=1,NC
422 TT(J)=R(I,J)
    DO 424 K=1,NC
    RT=0.
    DO 430 L=1,NC
430 RT=RT + TT(L)*VEC(L,K)
424 R(I,K)=RT
420 CONTINUE
C
C    CALCULATE ERROR BOUNDS AND CHECK FOR CONVERGENCE OF EIGENVALUES
C
    DO 380 I=1,NC
    VDOT=0.
    DO 382 J=1,NC
382 VDOT=VDOT + VEC(J.I)*VEC(J.I)
    EIGV2=EIGV(I)*EIGV(I)
    DIF=VDOT - EIGV2
    RDIF=MAX(DIF,TOLJ2*EIGV2)/VDOT
    RDIF=SQRT(RDIF)
    RTOLV(I)=RDIF
380 CONTINUE
    IF (IFPR.EQ.0 .AND. ICONV.EQ.0) GO TO 385
    WRITE (IOUT,1050)
    WRITE (IOUT,1005) (RTOLV(I),I=1,NC)
385 IF (ICONV.GT.0) GO TO 500
C
    DO 390 I=1,NROOT
    IF (RTOLV(I).GT.RTOL) GO TO 400
390 CONTINUE
    WRITE (IOUT,1060) RTOL
    ICONV=1
    GO TO 100
400 IF (NITE.LT.NITEM) GO TO 100
    WRITE (IOUT,1070)
    ICONV=2
    IFSS=0
    GO TO 100
C
C -- - END OF ITERATION LOOP
C
500 WRITE (IOUT,1100)
    WRITE (IOUT,1006) (EIGV(I),I=1,NROOT)
    WRITE (IOUT,1110)
    DO 530 J=1,NROOT
530 WRITE (IOUT,1005) (R(K,J),K=1,NN) 
```

<!-- source-page: 987 -->

```csv
C
C CALCULATE AND PRINT ERROR MEASURES
C
REWIND NSTIF
READ (NSTIF) A
C
DO 580 L=1,NROOT
RT-EIGV(L)
CALL MULT(TT,A,R(1,L),MAXA,NN,NWK)
VNORM=0.
DO 590 I=1,NN
590 VNORM=VNORM + TT(I)*TT(I)
CALL MULT(W,B,R(1,L),MAXA,NN,NWM)
WNORM=0.
DO 600 I=1,NN
TT(I)=TT(I) - RT*W(I)
600 WNORM=WNORM + TT(I)*TT(I)
VNORM=SQRT(VNORM)
WNORM=SQRT(WNORM)
D(L)=WNORM/VNORM
580 CONTINUE
WRITE (IOUT,1115)
WRITE (IOUT,1005) (D(I),I=1,NROOT)
C
APPLY STURM SEQUENCE CHECK
C
IF (IFSS.EQ.0) GO TO 900
CALL SCHECK (EIGV,RTOLV,BUP,BLO,BUPC,D,NC,NEI,RTOL,SHIFT,IOUT)
C
WRITE (IOUT,1120) SHIFT
C
SHIFT MATRIX A
C
REWIND NSTIF
READ (NSTIF) A
IF (NWM.GT.NN) GO TO 645
DO 640 I=1,NN
II=MAXA(I)
640 A(II)=A(II) - B(I)*SHIFT
GO TO 660
645 DO 650 I=1,NWK
650 A(I)=A(I) - B(I)*SHIFT
C
FACTORIZE SHIFTED MATRIX
C
660 ISH=1
CALL DECOMP (A,MAXA,NN,ISH,IOUT)
C
COUNT NUMBER OF NEGATIVE DIAGONAL ELEMENTS
C
NSCH=0
DO 664 I=1,NN
II=MAXA(I)
IF (A(II).LT.0.) NSCH=NSCH + 1
664 CONTINUE
IF (NSCH.EQ.NEI) GO TO 670
NMIS=NSCH - NEI
WRITE (IOUT,1130) NMIS
GO TO 900
670 WRITE (IOUT,1140) NSCH
GO TO 900
C
800 STOP
900 RETURN
C
1002 FORMAT ( ' ',10F10.0)
1005 FORMAT ( ' ',12E11.4)
1006 FORMAT ( ' ',6E22.14)
1007 FORMAT (///,' STOP, NC IS LARGER THAN THE NUMBER OF MASS ',1 'DEGREES OF FREEDOM')
SSP00282
SSP00283
SSP00284
SSP00285
SSP00286
SSP00287
SSP00288
SSP00289
SSP00290
SSP00291
SSP00292
SSP00293
SSP00294
SSP00295
SSP00296
SSP00297
SSP00298
SSP00299
SSP00300
SSP00301
SSP00302
SSP00303
SSP00304
SSP00305
SSP00306
SSP00307
SSP00308
SSP00309
SSP00310
SSP00311
SSP00312
SSP00313
SSP00314
SSP00315
SSP00316
SSP00317
SSP00318
SSP00319
SSP00320
SSP00321
SSP00322
SSP00323
SSP00324
SSP00325
SSP00326
SSP00327
SSP00328
SSP00329
SSP00330
SSP00331
SSP00332
SSP00333
SSP00334
SSP00335
SSP00336
SSP00337
SSP00338
SSP00339
SSP00340
SSP00341
SSP00342
SSP00343
SSP00344
SSP00345
SSP00346
SSP00347
SSP00348
SSP00349
SSP00350
SSP00351 
```

<!-- source-page: 988 -->

```csv
1008 FORMAT (///,' DEGREES OF FREEDOM EXCITED BY UNIT STARTING ', SSP00352
1 'ITERATION VECTORS') SSP00353
1010 FORMAT (///,' I T E R A T I O N N U M B E R ',I8) SSP00354
1020 FORMAT (///,' PROJECTION OF A (MATRIX AR)') SSP00355
1030 FORMAT (///,' PROJECTION OF B (MATRIX BR)') SSP00356
1035 FORMAT (///,' EIGENVALUES OF AR-LAMBDA*BR') SSP00357
1040 FORMAT (///,' AR AND BR AFTER JACOBI DIAGONALIZATION') SSP00358
1050 FORMAT (///,' ERROR BOUNDS REACHED ON EIGENVALUES') SSP00359
1060 FORMAT (///,' CONVERGENCE REACHED FOR RTOL ',E10.4) SSP00360
1070 FORMAT (' *** NO CONVERGENCE IN MAXIMUM NUMBER OF ITERATIONS', SSP00361
1 ' PERMITTED',/,
2 ' WE ACCEPT CURRENT ITERATION VALUES',/,
3 ' THE STURM SEQUENCE CHECK IS NOT PERFORMED') SSP00364
1100 FORMAT (///,' THE CALCULATED EIGENVALUES ARE') SSP00365
1115 FORMAT (///,' ERROR MEASURES ON THE EIGENVALUES') SSP00366
1110 FORMAT (///,' THE CALCULATED EIGENVECTORS ARE',/)
1120 FORMAT (///,' CHECK APPLIED AT SHIFT ',E22.14) SSP00368
1130 FORMAT (///,' THERE ARE ',I8,' EIGENVALUES MISSING') SSP00369
1140 FORMAT (///,' WE FOUND THE LOWEST ',I8,' EIGENVALUES') SSP00370
C
END
SUBROUTINE DECOMP (A,MAXA,NN,ISH,IOUT) SSP00372
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
P R O G R A M
TO CALCULATE (L)*(D)*(L)(T) FACTORIZATION OF
STIFFNESS MATRIX
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SSP00381
DIMENSION A(1),MAXA(1) SSP00382
IF (NN.EQ.1) GO TO 900 SSP00383
C
DO 200 N=1,NN
KN=MAXA(N) SSP00384
KL=KN + 1 SSP00385
KU=MAXA(N+1) - 1 SSP00386
KH=KU - KL SSP00387
IF (KH) 304,240,210 SSP00388
210 K=N - KH SSP00389
IC=0 SSP00390
KLT=KU SSP00391
DO 260 J=1,KH SSP00391
IC=IC + 1 SSP00392
KLT=KLT - 1 SSP00393
KI=MAXA(K) SSP00394
ND=MAXA(K+1) - KI - 1 SSP00395
IF (ND) 260,260,270 SSP00396
270 KK=MIN0(IC,ND) SSP00397
C=0. SSP00398
DO 280 L=1,KK SSP00399
280 C=C + A(KI+L)*A(KLT+L) SSP00400
A(KLT)=A(KLT) - C SSP00401
260 K=K + 1 SSP00402
240 K=N SSP00403
B=0. SSP00404
DO 300 KK=KL,KU SSP00405
K=K - 1 SSP00406
KI=MAXA(K) SSP00407
C=A(KK)/A(KI) SSP00408
IF (ABS(C).LT.1.E07) GO TO 290 SSP00409
WRITE (IOUT,2010) N,C SSP00410
GO TO 800 SSP00411
290 B=B + C*A(KK) SSP00412
300 A(KK)=C SSP00413
A(KN)=A(KN) - B SSP00414
304 IF (A(KN)) 310,310,200 SSP00415
310 IF (ISH.EQ.0) GO TO 320 SSP00416
IF (A(KN).EQ.0.) A(KN)=-1.E-16 SSP00417
300 K=MIN0(IC,ND) SSP00418
300 K=MIN0(IC,ND) SSP00419
300 K=MIN0(IC,ND) SSP00420
300 K=MIN0(IC,ND) SSP00421 
```

<!-- source-page: 989 -->

```txt
GO TO 200
320 WRITE (IOUT,2000) N,A(KN)
GO TO 800
200 CONTINUE
GO TO 900
C
800 STOP
900 RETURN
C
2000 FORMAT (/// STOP - STIFFNESS MATRIX NOT POSITIVE DEFINITE',///,
1    'NONPOSITIVE PIVOT FOR EQUATION ',I8,//,
2    'PIVOT = ',E20.12)
2010 FORMAT (/// STOP - STURM SEQUENCE CHECK FAILED BECAUSE OF',
1    'MULTIPLIER GROWTH FOR COLUMN NUMBER ',I8,//,
2    'MULTIPLIER = ',E20.8)
END
SUBROUTINE REDBAK (A,V,MAXA,NN)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(1),V(1),MAXA(1)
C
DO 400 N=1,NN
KL=MAXA(N) + 1
KU=MAXA(N+1) - 1
IF (KU-KL) 400,410,410
410 K=N
C=0.
DO 420 KK=KL,KU
K=K - 1
420 C=C + A(KK)*V(K)
V(N)=V(N) - C
400 CONTINUE
C
DO 480 N=1,NN
K=MAXA(N)
480 V(N)=V(N)/A(K)
IF (NN.EQ.1) GO TO 900
N=NN
DO 500 L=2,NN
KL=MAXA(N) + 1
KU=MAXA(N+1) - 1
IF (KU-KL) 500,510,510
510 K=N
DO 520 KK=KL,KU
K=K - 1
520 V(K)=V(K) - A(KK)*V(N)
500 N=N - 1
C
900 RETURN
END
SUBROUTINE MULT (TT,B,RR,MAXA,NN,NWM)
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
320
32
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
3
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
T
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
P
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
S
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
M
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
G
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
N
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
O
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
D
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
U
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
I
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
IC
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
R
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
B
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
E
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
6
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
0
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
L
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
F
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
Cc
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
V
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
CI
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
CC
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
NC
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
CP
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
CT
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
H
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
J
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
No
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
<|content_end|> 
```

<!-- source-page: 990 -->

```csv
GO TO 900
SSP00492
SSP00493
20 DO 40 I=1,NN
SSP00494
40 TT(I)=0.
SSP00495
DO 100 I=1,NN
SSP00496
KL=MAXA(I)
SSP00497
KU=MAXA(I+1) - 1
SSP00498
II=I + 1
SSP00499
CC=RR(I)
SSP00500
DO 100 KK=KL,KU
SSP00501
II=II - 1
SSP00502
100 TT(II)=TT(II) + B(KK)*CC
SSP00503
IF (NN.EQ.1) GO TO 900
SSP00504
DO 200 I=2,NN
SSP00505
KL=MAXA(I) + 1
SSP00506
KU=MAXA(I+1) - 1
SSP00507
IF (KU-KL) 200,210,210
SSP00508
210 II=I
SSP00509
AA=0.
SSP00510
DO 220 KK=KL,KU
SSP00511
II=II - 1
SSP00512
220 AA=AA + B(KK)*RR(II)
SSP00513
TT(I)=TT(I) + AA
SSP00514
200 CONTINUE
SSP00515
900 RETURN
SSP00516
END
SSP00517
SUBROUTINE SCHECK (EIGV,RTOLV,BUP,BLO,BUPC,NEIV,NC,NEI,RTOL,
1 SHIFT,IOUT)
SSP00518
SSP00519
SSP00520
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
P R O G R A M
SSP00523
C . TO EVALUATE SHIFT FOR STURM SEQUENCE CHECK
SSP00524
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SSP00527
DIMENSION EIGV(NC),RTOLV(NC),BUP(NC),BLO(NC),BUPC(NC),NEIV(NC)
SSP00528
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION EIGV(NC),RTOLV(NC),BUP(NC),BLO(NC),BUPC(NC),NEIV(NC)
SSP00529
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION EIGV(NC),RTOLV(NC),BUP(NC),BLO(NC),BUPC(NC),NEIV(NC)
SSP00530
FTOL=0.01
SSP00531
C DO 100 I=1,NC
SSP00532
BUP(I)=EIGV(I)*(1.+FTOL)
SSP00533
100 BLO(I)=EIGV(I)*(1.-FTOL)
SSP00534
NROOT=0
SSP00535
DO 120 I=1,NC
SSP00536
120 IF (RTOLV(I).LT.RTOL) NROOT=NROOT + 1
SSP00537
IF (NROOT.GE.1) GO TO 200
SSP00538
WRITE (IOUT,1010)
SSP00539
GO TO 800
SSP00540
SSP00541
C FIND UPPER BOUNDS ON EIGENVALUE CLUSTERS
SSP00542
C FIND UPPER BOUNDS ON EIGENVALUE CLUSTERS
SSP00543
200 DO 240 I=1,NROOT
SSP00544
240 NEIV(I)=1
SSP00545
IF (NROOT.NE.1) GO TO 260
SSP00546
BUPC(1)=BUP(1)
SSP00547
LM=1
SSP00548
L=1
SSP00549
I=2
SSP00550
GO TO 295
SSP00551
260 L=1
SSP00552
I=2
SSP00553
270 IF (BUP(I-1).LE.BLO(I)) GO TO 280
SSP00554
NEIV(L)=NEIV(L) + 1
SSP00555
I=I + 1
SSP00556
IF (I.LE.NROOT) GO TO 270
SSP00557
280 BUPC(L)=BUP(I-1)
IF (I.GT.NROOT) GO TO 290
SSP00558
L=L + 1
SSP00559
SSP00560
SSP00561 
```