-
Notifications
You must be signed in to change notification settings - Fork 17
/
newnumal5p2.txt
21227 lines (16951 loc) · 773 KB
/
newnumal5p2.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 1
AUTHORS: T.J.DEKKER AND W.HOFFMANN.
CONTRIBUTORS: W.HOFFMANN, J.G.VERWER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 730716.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS FOUR PROCEDURES FOR CALCULATING EIGENVALUES
OR EIGENVECTORS OF A SYMMETRIC TRIDIAGONAL MATRIX.
VALSYMTRI CALCULATES ALL, OR SOME CONSECUTIVE, EIGENVALUES OF A
SYMMETRIC TRIDIAGONAL MATRIX BY MEANS OF LINEAR INTERPOLATION USING
A STURM SEQUENCE;
VECSYMTRI CALCULATES THE CORRESPONDING EIGENVECTORS BY MEANS OF
INVERSE ITERATION.
QRIVALSYMTRI CALCULATES ALL EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
MATRIX BY MEANS OF QR ITERATION;
QRISYMTRI CALCULATES THE EIGENVECTORS AS WELL.
WHEN ALL EIGENVALUES HAVE TO BE CALCULATED, QRIVALSYMTRI IS
PREFERABLE WITH RESPECT TO THE RUNNING TIME; WHEN THE EIGENVECTORS
ALSO HAVE TO BE CALCULATED, INVERSE ITERATION IS PREFERABLE.
KEYWORDS:
EIGENVALUES,
EIGENVECTORS,
TRIDIAGONAL MATRIX,
STURM-SEQUENCE,
INVERSE ITERATION,
QR ITERATION.
1SECTION 3.3.1.1.1 (DECEMBER 1975) PAGE 2
SUBSECTION: VALSYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM);
"VALUE" N, N1, N2; "INTEGER" N, N1, N2;
"ARRAY" D, BB, VAL, EM;
"CODE" 34151;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N];
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
BB: <ARRAY IDENTIFIER>;
"ARRAY" BB[1:N-1];
ENTRY: THE SQUARES OF THE CODIAGONAL ELEMENTS OF THE
SYMMETRIC TRIDIAGONAL MATRIX;
N1,N2: <ARITHMETIC EXPRESSION>;
THE SERIAL NUMBER OF THE FIRST AND LAST EIGENVALUE TO BE
CALCULATED, RESPECTIVELY;
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[N1:N2];
EXIT: THE N2-N1+1 CALCULATED CONSECUTIVE EIGENVALUES IN
NONINCREASING ORDER;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:3];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[1], AN UPPERBOUND FOR THE MODULI OF THE
EIGENVALUES OF THE GIVEN MATRIX,
EM[2], A RELATIVE TOLERANCE FOR THE EIGENVALUES;
EXIT: EM[3], THE TOTAL NUMBER OF ITERATIONS USED FOR
CALCULATING THE EIGENVALUES.
PROCEDURES USED:
ZEROIN = CP34150.
RUNNING TIME:
DEPENDS STRONGLY ON THE DISTANCE OF SUCCESSIVE EIGENVALUES.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
LET T DENOTE THE GIVEN SYMMETRIC TRIDIAGONAL MATRIX OF ORDER N AND
I THE IDENTITY MATRIX. THE EIGENVALUES OF T ARE THE ZEROES OF THE
N-TH DEGREE POLYNOMIAL P(N,X) = DET(T - X*I). INSTEAD OF SEARCHING
FOR THE ZEROES OF P(N,X) WE LOOK FOR THE ZEROES OF THE FUNCTION
F(N,X) = P(N,X) / P(N-1,X). MAINTAINING A LOWER BOUND FOR
ABS(P(N-1,X)) WE DO AVOID OVERFLOW OF THE REAL NUMBER CAPACITY IN
THE COMPUTATION OF F(N,X). THIS FUNCTION CAN BE CALCULATED AS
FOLLOWS:
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 3
F(1,X) = D[1] - X,
F(K,X) = D[K] - X - BB[K-1] /
("IF" ABS(F(K-1,X)) > MACHTOL "THEN" F(K-1,X)
"ELSE" "IF" F(K-1,X) <= 0 "THEN" -MACHTOL
"ELSE" MACHTOL), K = 2, . . . ,N,
WHERE MACHTOL EQUALS EM[0] * EM[1].
USING THE STURM SEQUENCE PROPERTY OF (F(K,X)), K=1,2,...,N, WE CAN
LOCATE THE DESIRED EIGENVALUES BY MEANS OF THE PROCEDURE ZEROIN
(SECTION 5.1.1.1). FOR FURTHER DETAILS SEE REF[1], REF[2].
SUBSECTION: VECSYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" VECSYMTRI(D, B, N, N1, N2, VAL, VEC, EM);
"VALUE" N, N1, N2; "INTEGER" N, N1, N2;
"ARRAY" D, B, VAL, VEC, EM;
"CODE" 34152;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N],
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
B: <ARRAY IDENTIFIER>;
"ARRAY" B[1:N];
ENTRY: THE CODIAGONAL OF THE SYMMETRIC TRIDIAGONAL MATRIX
FOLLOWED BY AN ADDITIONAL ELEMENT 0;
N1, N2: <ARITHMETIC EXPRESSION>;
LOWER AND UPPER BOUND OF THE ARRAY VAL (SEE ALSO METHOD AND
PERFORMANCE);
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[N1:N2];
ENTRY: A ROW OF NONINCREASING EIGENVALUES AS DELIVERED BY
VALSYMTRI;
VEC: <ARRAY IDENTIFIER>;
"ARRAY" VEC[1:N,N1:N2];
EXIT: THE EIGENVECTORS CORRESPONDING WITH THE GIVEN
EIGENVALUES (SEE ALSO METHOD AND PERFORMANCE);
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:9];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[1], A NORM OF THE GIVEN MATRIX,
EM[4], THE ORTHOGONALISATION PARAMETER (SEE ALSO
METHOD AND PERFORMANCE),
EM[6], THE RELATIVE TOLERANCE FOR THE EIGENVECTORS,
EM[8], THE MAXIMUM NUMBER OF ITERATIONS ALLOWED
FOR THE CALCULATION OF EACH EIGENVECTOR;
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 4
EXIT: EM[5], THE NUMBER OF EIGENVECTORS INVOLVED IN THE
LAST GRAM-SCHMIDT ORTHOGONALISATION (SEE
METHOD AND PERFORMANCE),
EM[7], THE MAXIMUM EUCLIDEAN NORM OF THE RESIDUES,
EM[9], THE LARGEST NUMBER OF ITERATIONS PERFORMED
FOR THE CALCULATION OF SOME EIGENVECTOR (SEE
METHOD AND PERFORMANCE).
PROCEDURES USED:
VECVEC = CP34010,
TAMVEC = CP34012,
ELMVECCOL= CP34021.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: FIVE AUXILIARY ONE-DIMENSIONAL REAL ARRAYS
AND ONE BOOLEAN ARRAY, ALL OF LENGTH N, ARE USED.
RUNNING TIME: THE PROCESS IS OF ORDER N FOR EACH EIGENVECTOR.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
AN EIGENVECTOR OF A SYMMETRIC TRIDIAGONAL MATRIX T, CORRESPONDING
TO AN EIGENVALUE LAMBDA, IS CALCULATED BY MEANS OF INVERSE
ITERATION; I.E. STARTING FROM SOME INITIAL VECTOR X, THE LINEAR
SYSTEM (T - LAMBDA * I)Y = X IS SOLVED ITERATIVELY, THE SOLUTION Y,
DIVIDED BY ITS EUCLIDEAN NORM REPLACING X EACH TIME.
IF THE DISTANCE BETWEEN SOME APPROXIMATE EIGENVALUES IS SMALLER
THAN MACHTOL (=EM[0] * EM[1]), THEN THEY ARE SLIGHTLY MODIFIED SUCH
THAT THE DISTANCE BETWEEN THEM EQUALS MACHTOL. IF THE DISTANCE
BETWEEN SOME EIGENVALUES IS SMALLER THAN THE ORTHOGONALISATION
PARAMETER (=EM[4]) TIMES EM[1], THEN IN EACH ITERATION STEP GRAM-
SCHMIDT ORTHOGONALISATION IS CARRIED OUT, SO THAT THE EIGENVECTORS
OBTAINED ARE ORTHOGONAL WITHIN WORKING PRECISION. THE ITERATION
ENDS AS SOON AS EITHER THE EUCLIDEAN NORM OF THE RESIDUE IS SMALLER
THAN EM[1] * EM[6], OR THE MAXIMUM ALLOWED NUMBER OF ITERATIONS
(=EM[8]) HAS BEEN PERFORMED. IN THE LATTER CASE EM[9]:= EM[8] + 1.
IF N1 > 1, THEN VECSYMTRI SHOULD BE PRECEDED BY ONE OR MORE CALLS
OF VECSYMTRI PRODUCING A NUMBER OF EIGENVECTORS CORRESPONDING TO
THE PRECEDING EIGENVALUES. MOREOVER ONE MUST GIVE EM[5], AS
PRODUCED BY THE LAST CALL OF VECSYMTRI; THE K-TH TO N2-TH
EIGENVALUES, WHERE K = N1 - EM[5], MUST BE GIVEN IN ARRAY VAL[K:N2]
IN MONOTONICALLY NONINCREASING ORDER (THE K-TH TO (N1-1)-TH
EIGENVALUES BEING NEEDED FOR THE MODIFYING MENTIONED ABOVE), AND
THE CORRESPONDING EIGENVECTORS UP TO THE (N1-1)-TH (WHICH ARE
NEEDED FOR THE GRAM-SCHMIDT ORTHOGONALISATION) IN THE CORRESPONDING
COLUMNS OF ARRAY VEC[1:N,K:N2].
THE TOLERANCES SHOULD SATISFY: EM[0] ( < EM[2]) < EM[6] AND
EM[4] >= EM[0] / EM[6]. FOR FURTHER DETAILS SEE REF[1].
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 5
SUBSECTION: QRIVALSYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"INTEGER" "PROCEDURE" QRIVALSYMTRI(D, BB, N, EM);
"VALUE" N; "INTEGER" N; "ARRAY" D, BB, EM;
"CODE" 34160;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N];
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
EXIT: THE EIGENVALUES OF THE MATRIX IN SOME ARBITRARY
ORDER;
BB: <ARRAY IDENTIFIER>;
"ARRAY" BB[1:N];
ENTRY: THE SQUARES OF THE CODIAGONAL ELEMENTS OF THE
SYMMETRIC TRIDIAGONAL MATRIX FOLLOWED BY AN
ADDITONAL ELEMENT O;
EXIT: THE SQUARES OF THE CODIAGONAL ELEMENTS OF THE
SYMMETRIC TRIDIAGONAL MATRIX RESULTING FROM THE QR
ITERATION;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:5];
ENTRY: EM[0], THE MACHINE PRECISION;
EM[1], A NORM OF THE GIVEN MATRIX;
EM[2], A RELATIVE TOLERANCE FOR THE EIGENVALUES;
EM[4], THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;
EXIT: EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE CODIAGONAL
ELEMENTS NEGLECTED;
EM[5], THE NUMBER OF ITERATIONS PERFORMED.
MOREOVER:
QRIVALSYMTRI:= THE NUMBER OF EIGENVALUES NOT CALCULATED.
PROCEDURES USED: NONE.
RUNNING TIME: THE PROCESS IS OF ORDER N SQUARED.
LANGUAGE: ALGOL 60.
1SECTION 3.3.1.1.1 (DECEMBER 1975) PAGE 6
METHOD AND PERFORMANCE:
IN QRIVALSYMTRI THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX
ARE CALCULATED BY MEANS OF QR-ITERATION. FOR THIS PROCEDURE WE USED
ESSENTIALLY THE SQUARE-ROOT-FREE VERSION OF THE QR ALGORITHM DUE TO
REINSCH[3].
IN ADDITION TO THE RELATIVE ERROR, WHICH IS SUPPOSED TO BE BOUNDED
BY EM[1] * EM[2] (I.E. MATRIX NORM TIMES RELATIVE TOLERANCE), THE
CALCULATED EIGENVALUES HAVE AN ABSOLUTE ERROR WHICH IS BOUNDED BY
BY EM[0] * EM[1] (I.E. MACHINE PRECISION TIMES MATRIX NORM).
IN PARTICULAR, WHEN SOME EIGENVALUES ARE VERY SMALL COMPARED TO THE
MATRIX NORM, THE ACCURACY OF THE CALCULATED EIGENVALUES CAN BE
INCREASED BY GIVING EM[0] A (POSITIVE) VALUE WHICH IS LESS THAN
THE MACHINE PRECISION.
A PARTICULAR CHOICE OF EM[0] IS HARMLESS FOR THE PROCEDURE
PROVIDED THAT FOR EACH I THE CALCULATION OF BB[I] / EM[0] ** 2
CAUSES NO OVERFLOW AND THE CALCULATION OF (EM[0] * EM[1]) ** 2
CAUSES NO UNDERFLOW.
ONE SHOULD NOTICE THAT THE NUMBER OF QR ITERATIONS INCREASES BY A
SMALLER CHOICE OF EM[0].
FOR FURTHER DETAILS SEE [2], [3].
SUBSECTION: QRISYMTRI.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"INTEGER" "PROCEDURE" QRISYMTRI(A, N, D, B, BB, EM);
"VALUE" N; "INTEGER" N; "ARRAY" A, D, B, BB, EM;
"CODE" 34161;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1:N];
ENTRY: THE MAIN DIAGONAL OF THE SYMMETRIC TRIDIAGONAL
MATRIX;
EXIT: THE EIGENVALUES OF THE MATRIX IN SOME ARBITRARY
ORDER;
B: <ARRAY IDENTIFIER>;
"ARRAY" B[1:N];
ENTRY: THE CODIAGONAL OF THE SYMMETRIC TRIDIAGONAL MATRIX
FOLLOWED BY AN ADDITIONAL ELEMENT 0;
EXIT: THE CODIAGONAL OF THE SYMMETRIC TRIDIAGONAL MATRIX
RESULTING FROM THE QR ITERATION, FOLLOWED BY AN
ADDITIONAL ELEMENT 0;
BB: <ARRAY IDENTIFIER>;
"ARRAY" BB[1:N];
ENTRY: THE SQUARED CODIAGONAL ELEMENTS OF THE SYMMETRIC
TRIDIAGONAL MATRIX, FOLLOWED BY AN ADDITIONAL
ELEMENT O;
EXIT: THE SQUARED CODIAGONAL ELEMENTS OF THE SYMMETRIC
TRIDIAGONAL MATRIX RESULTING FROM THE QR ITERATION;
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 7
A: <ARRAY IDENTIFIER>;
"ARRAY" A[1:N,1:N];
ENTRY: SOME MATRIX S, SAY, (POSSIBLY THE IDENTITY MATRIX);
EXIT: THE EIGENVECTORS OF THE ORIGINAL SYMMETRIC
TRIDIAGONAL MATRIX, PREMULTIPLIED BY S (SEE METHOD
AND PERFORMANCE);
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:5];
ENTRY: EM[0], THE MACHINE PRECISION;
EM[1], A NORM OF THE GIVEN MATRIX;
EM[2], A RELATIVE TOLERANCE FOR THE QR ITERATION;
EM[4], THE MAXIMUM ALLOWED NUMBER OF ITERATIONS;
EXIT: EM[3], THE MAXIMUM ABSOLUTE VALUE OF THE CODIAGONAL
ELEMENTS NEGLECTED;
EM[5], THE NUMBER OF ITERATIONS PERFORMED.
MOREOVER:
QRISYMTRI:= THE NUMBER OF EIGENVALUES AND -VECTORS NOT
CALCULATED.
PROCEDURES USED:
ROTCOL = CP34040.
RUNNING TIME: THE PROCESS IS OF ORDER N CUBED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
IN QRISYMTRI THE EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC
TRIDIAGONAL MATRIX ARE COMPUTED SIMULTANEOUSLY.
IN MOST APPLICATIONS QRISYMTRI IS USED IN THE COMPUTATION OF
EIGENVALUES AND -VECTORS OF A GENERAL SYMMETRIC MATRIX (SEE QRISYM
SECTION 3.3.1.1.2); IN THAT CASE ARRAY A IS INITIALLY GIVEN THE
VALUE OF THE TRANSFORMING MATRIX (TFMPREVEC SECTION 3.2.1.2.1.1).
FOR THE COMPUTATION OF EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC
TRIDIAGONAL MATRIX, ARRAY A HAS TO BE INITIALIZED TO THE IDENTITY
MATRIX. THE AVERAGE NUMBER OF ITERATIONS IS ABOUT 3N. WHEN THE
PROCESS IS COMPLETED WITHIN EM[4] ITERATIONS, THEN QRISYMTRI:= 0;
OTHERWISE QRISYMTRI:= THE NUMBER, K, OF EIGENVALUES NOT CALCULATED,
EM[5]:= EM[4] + 1 AND ONLY THE LAST N - K ELEMENTS OF D AND THE
LAST N - K COLUMNS OF A ARE APPROXIMATE EIGENVALUES AND -VECTORS
RESPECTIVELY, OF THE ORIGINAL MATRIX.
FOR FURTHER DETAILS SEE REF[1], REF[2].
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 8
REFERENCES:
[1] DEKKER, T.J. AND HOFFMANN, W.
ALGOL 60 PROCEDURES IN NUMERICAL ALGEBRA, PART 2,
MATHEMATICAL CENTRE TRACTS 23,
MATHEMATISCH CENTRUM, AMSTERDAM, 1968;
[2] WILKINSON, J.H.
THE ALGEBRAIC EIGENVALUE PROBLEM,
CLARENDON PRESS, OXFORD 1965.
[3] REINSCH, CHR.H.
A STABLE, RATIONAL QR ALGORITHM FOR THE COMPUTATION OF THE
EIGENVALUES OF AN HERMITIAN, TRIDIAGONAL MATRIX.
MATH. OF COMP. VOL 25(1971) PP. 591-597.
EXAMPLE OF USE:
THE FIRST AND SECOND EIGENVALUE IN MONOTONICALLY NON-INCREASING
ORDER AND THE CORRESPONDING EIGENVECTORS OF T, WITH N = 4 AND
T[I,J] = "IF" I = J "THEN" 2 "ELSE" "IF" ABS(I - J) = 1 "THEN" - 1
"ELSE" 0, MAY BE OBTAINED BY THE FOLLOWING PROGRAM:
"BEGIN"
"INTEGER" J;
"ARRAY" B, D[1:4], BB[1:3], VAL[1:2], EM[0:9], VEC[1:4,1:2];
EM[0]:= "-14; EM[1]:= 4; EM[2]:= "-12;
EM[4]:= "-3; EM[6]:= "-10; EM[8]:= 5;
"FOR" J:= 1, 2, 3, 4 "DO" D[J]:= 2; B[4]:= 0;
"FOR" J:= 1, 2, 3 "DO"
"BEGIN" BB[J]:= 1; B[J]:= -1 "END";
VALSYMTRI(D, BB, 4, 1, 2, VAL, EM);
VECSYMTRI(D, B, 4, 1, 2, VAL, VEC, EM);
OUTPUT(61, "("2(+.13D"+2D, 2B), 2/")", VAL[1], VAL[2]);
"FOR" J:= 1, 2, 3, 4 "DO"
OUTPUT(61, "("2(+.13D"+2D, 2B), /")", VEC[J,1], VEC[J,2]);
OUTPUT(61, "("/, .2D"+2D, /, 3(2ZD, /)")",
EM[7], EM[3], EM[5], EM[9])
"END"
THE PROGRAM DELIVERS:
THE EIGENVALUES: +.3618033988751"+01 +.2618033988750"+01
THE EIGENVECTORS: +.3717480344602"+00 +.6015009550075"+00
+.6015009550075"+00 +.3717480344602"+00
+.6015009550075"+00 -.3717480344602"+00
+.3717480344602"+00 -.6015009550075"+00
EM[7] = .15"-11
EM[3] = 24
EM[5] = 1
EM[9] = 1 .
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 9
SOURCE TEXT(S):
0"CODE" 34151;
"COMMENT" MCA 2311;
"PROCEDURE" VALSYMTRI(D, BB, N, N1, N2, VAL, EM);
"VALUE" N, N1, N2;
"INTEGER" N, N1, N2; "ARRAY" D, BB, VAL, EM;
"BEGIN" "INTEGER" K, COUNT;
"REAL" MAX, X, Y, MACHEPS, NORM, RE, MACHTOL, UB, LB, LAMBDA;
"REAL" "PROCEDURE" STURM;
"BEGIN" "INTEGER" P, I; "REAL" F;
COUNT:= COUNT + 1;
P:= K; F:= D[1] - X;
"FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"
"BEGIN" "IF" F <= 0 "THEN"
"BEGIN" P:= P + 1;
"IF" P > N "THEN" "GOTO" OUT
"END"
"ELSE" "IF" P < I - 1 "THEN"
"BEGIN" LB:= X; "GOTO" OUT "END";
"IF" ABS(F) < MACHTOL "THEN"
F:= "IF" F <= 0 "THEN" - MACHTOL "ELSE" MACHTOL;
F:= D[I] - X - BB[I - 1] / F
"END";
"IF" P = N "OR" F <= 0 "THEN"
"BEGIN" "IF" X < UB "THEN" UB:= X "END" "ELSE" LB:= X;
OUT: STURM:= "IF" P = N "THEN" F "ELSE" (N - P) * MAX
"END" STURM;
MACHEPS:= EM[0]; NORM:= EM[1]; RE:= EM[2];
MACHTOL:= NORM * MACHEPS; MAX:= NORM / MACHEPS; COUNT:= 0;
UB:= 1.1 * NORM; LB:= - UB; LAMBDA:= UB;
"FOR" K:= N1 "STEP" 1 "UNTIL" N2 "DO"
"BEGIN" X:= LB; Y:= UB; LB:= -1.1 * NORM;
ZEROIN(X, Y, STURM, ABS(X) * RE + MACHTOL);
VAL[K]:= LAMBDA:= "IF" X > LAMBDA "THEN" LAMBDA "ELSE" X;
"IF" UB > X "THEN" UB:= "IF" X > Y "THEN" X "ELSE" Y
"END";
EM[3]:= COUNT
"END" VALSYMTRI
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 10
;
"EOP"
0"CODE" 34152;
"COMMENT" MCA 2312;
"PROCEDURE" VECSYMTRI(D, B, N, N1, N2, VAL, VEC, EM);
"VALUE" N, N1, N2;
"INTEGER" N, N1, N2; "ARRAY" D, B, VAL, VEC, EM;
"BEGIN" "INTEGER" I, J, K, COUNT, MAXCOUNT, COUNTLIM, ORTH, IND;
"REAL" BI, BI1, U, W, Y, MI1, LAMBDA, OLDLAMBDA, ORTHEPS,
VALSPREAD, SPR, RES, MAXRES, OLDRES, NORM, NEWNORM, OLDNORM,
MACHTOL, VECTOL;
"ARRAY" M, P, Q, R, X[1:N];
"BOOLEAN" "ARRAY" INT[1:N];
NORM:= EM[1]; MACHTOL:= EM[0] * NORM; VALSPREAD:= EM[4] * NORM;
VECTOL:= EM[6] * NORM; COUNTLIM:= EM[8]; ORTHEPS:= SQRT(EM[0]);
MAXCOUNT:= IND:= 0; MAXRES:= 0;
"IF" N1 > 1 "THEN"
"BEGIN" ORTH:= EM[5]; OLDLAMBDA:= VAL[N1 - ORTH];
"FOR" K:= N1 - ORTH + 1 "STEP" 1 "UNTIL" N1 - 1 "DO"
"BEGIN" LAMBDA:= VAL[K]; SPR:= OLDLAMBDA - LAMBDA;
"IF" SPR < MACHTOL "THEN" LAMBDA:= OLDLAMBDA - MACHTOL;
OLDLAMBDA:= LAMBDA
"END"
"END" "ELSE" ORTH:= 1;
"FOR" K:= N1 "STEP" 1 "UNTIL" N2 "DO"
"BEGIN" LAMBDA:= VAL[K]; "IF" K > 1 "THEN"
"BEGIN" SPR:= OLDLAMBDA - LAMBDA;
"IF" SPR < VALSPREAD "THEN"
"BEGIN" "IF" SPR < MACHTOL "THEN"
LAMBDA:= OLDLAMBDA - MACHTOL;
ORTH:= ORTH +1
"END" "ELSE" ORTH:= 1
"END";
COUNT:= 0; U:= D[1] - LAMBDA; BI:= W:= B[1];
"IF" ABS(BI) < MACHTOL "THEN" BI:= MACHTOL;
"FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
"BEGIN" BI1:= B[I + 1];
"IF" ABS(BI1) < MACHTOL "THEN" BI1:= MACHTOL;
"IF" ABS(BI) >= ABS(U) "THEN"
"BEGIN" MI1:= M[I + 1]:= U / BI; P[I]:= BI;
Y:= Q[I]:= D[I + 1] - LAMBDA; R[I]:= BI1;
U:= W - MI1 * Y; W:= - MI1 * BI1; INT[I]:= "TRUE"
"END"
"ELSE"
"BEGIN" MI1:= M[I + 1]:= BI / U; P[I]:= U; Q[I]:= W;
R[I]:= 0; U:= D[I + 1] - LAMBDA - MI1 * W;W:= BI1;
INT[I]:= "FALSE"
"END";
X[I]:= 1; BI:= BI1
"END" TRANSFORM
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 11
;
P[N]:= "IF" ABS(U) < MACHTOL "THEN" MACHTOL "ELSE" U;
Q[N]:= R[N]:= 0; X[N]:= 1; "GOTO" ENTRY;
ITERATE: W:= X[1];
"FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"
"BEGIN" "IF" INT[I - 1] "THEN"
"BEGIN" U:= W; W:= X[I - 1]:= X[I] "END"
"ELSE" U:= X[I]; W:= X[I]:= U - M[I] * W
"END" ALTERNATE;
ENTRY: U:= W:= 0;
"FOR" I:= N "STEP" -1 "UNTIL" 1 "DO"
"BEGIN" Y:= U; U:= X[I]:= (X[I] - Q[I] * U - R[I] * W) /
P[I]; W:= Y
"END" NEXT ITERATION;
NEWNORM:= SQRT(VECVEC(1, N, 0, X, X)); "IF" ORTH > 1"THEN"
"BEGIN" OLDNORM:= NEWNORM;
"FOR" J:= K - ORTH + 1 "STEP" 1 "UNTIL" K - 1 "DO"
ELMVECCOL(1, N, J, X, VEC, -TAMVEC(1, N, J, VEC, X));
NEWNORM:= SQRT(VECVEC(1, N, 0, X, X));
"IF" NEWNORM < ORTHEPS * OLDNORM "THEN"
"BEGIN" IND:= IND + 1; COUNT:= 1;
"FOR" I:= 1 "STEP" 1 "UNTIL" IND - 1,
IND + 1 "STEP" 1 "UNTIL" N "DO" X[I]:= 0;
X[IND]:= 1; "IF" IND = N "THEN" IND:= 0;
"GOTO" ITERATE
"END" NEW START
"END" ORTHOGONALISATION;
RES:= 1 / NEWNORM; "IF" RES > VECTOL "OR" COUNT = 0 "THEN"
"BEGIN" COUNT:= COUNT + 1; "IF" COUNT <= COUNTLIM "THEN"
"BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
X[I]:= X[I] * RES; "GOTO" ITERATE
"END"
"END";
"FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" VEC[I,K]:= X[I] * RES;
"IF" COUNT > MAXCOUNT "THEN" MAXCOUNT:= COUNT;
"IF" RES > MAXRES "THEN" MAXRES:= RES; OLDLAMBDA:= LAMBDA
"END";
EM[5]:= ORTH; EM[7]:= MAXRES; EM[9]:= MAXCOUNT
"END" VECSYMTRI
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 12
;
"EOP"
0"CODE" 34160;
"INTEGER" "PROCEDURE" QRIVALSYMTRI(D, BB, N, EM); "VALUE" N;
"INTEGER" N; "ARRAY" D, BB, EM;
"BEGIN" "INTEGER" I, I1, LOW, OLDLOW, N1, COUNT, MAX;
"REAL" BBTOL, BBMAX, BBI, BBN1, MACHTOL, DN, DELTA, F, NUM,
SHIFT, G, H, T, P, R, S, C, OLDG;
BBTOL:= (EM[2] * EM[1]) ** 2; MACHTOL:= EM[0] * EM[1];
MAX:= EM[4]; BBMAX:= 0; COUNT:= 0; OLDLOW:= N;
"FOR" N1:= N - 1 "WHILE" N > 0 "DO"
"BEGIN"
"FOR" I:= N, I - 1 "WHILE" ("IF" I >= 1 "THEN"
BB[I] > BBTOL "ELSE" "FALSE") "DO" LOW:= I;
"IF" LOW > 1 "THEN" "BEGIN" "IF" BB[LOW-1] > BBMAX "THEN"
BBMAX:= BB[LOW-1] "END";
"IF" LOW = N "THEN" N:= N1 "ELSE"
"BEGIN" DN:= D[N]; DELTA:= D[N1] - DN;
BBN1:= BB[N1];
"IF" ABS(DELTA) < MACHTOL "THEN" R:= SQRT(BBN1) "ELSE"
"BEGIN"
F:= 2 / DELTA; NUM:= BBN1 * F;
R:= -NUM / (SQRT(NUM * F + 1) + 1)
"END";
"IF" LOW = N1 "THEN"
"BEGIN" D[N]:= DN + R; D[N1]:= D[N1] - R; N:= N - 2
"END"
"ELSE"
"BEGIN" COUNT:= COUNT + 1;
"IF" COUNT > MAX "THEN" "GOTO" END;
"IF" LOW < OLDLOW "THEN"
"BEGIN" SHIFT:= 0; OLDLOW:= LOW "END"
"ELSE" SHIFT:= DN + R;
H:= D[LOW] - SHIFT;
"IF" ABS(H) < MACHTOL "THEN" H:= "IF" H <= 0 "THEN"
-MACHTOL "ELSE" MACHTOL;
G:= H; T:= G * H;
BBI:= BB[LOW]; P:= T + BBI; I1:= LOW;
"FOR" I:= LOW + 1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" S:= BBI / P; C:= T / P;
H:= D[I] - SHIFT - BBI / H;
"IF" ABS(H) < MACHTOL "THEN" H:= "IF" H <= 0
"THEN" -MACHTOL "ELSE" MACHTOL;
OLDG:= G; G:= H * C; T:= G * H;
D[I1]:= OLDG - G + D[I];
BBI:= "IF" I = N "THEN" 0 "ELSE" BB[I];
P:= T + BBI; BB[I1]:= S * P; I1:= I
"END";
D[N]:= G + SHIFT
"END" QRSTEP
"END"
"END";
END: EM[3]:= SQRT(BBMAX); EM[5]:= COUNT; QRIVALSYMTRI:= N
"END" QRIVALSYMTRI
1SECTION 3.3.1.1.1 (JULY 1974) PAGE 13
;
"EOP"
0"CODE" 34161;
"COMMENT" MCA 2321;
"INTEGER" "PROCEDURE" QRISYMTRI(A, N, D, B, BB, EM); "VALUE" N;
"INTEGER" N; "ARRAY" A, D, B, BB, EM;
"BEGIN" "INTEGER" I, J, J1, K, M, M1, COUNT, MAX;
"REAL" BBMAX, R, S, SIN, T, C, COS, OLDCOS, G, P, W, TOL, TOL2,
LAMBDA, DK1, A0, A1;
TOL:= EM[2] * EM[1]; TOL2:= TOL * TOL; COUNT:= 0; BBMAX:= 0;
MAX:= EM[4]; M:= N;
IN: K:= M; M1:= M - 1;
NEXT: K:= K - 1; "IF" K > 0 "THEN"
"BEGIN" "IF" BB[K] >= TOL2 "THEN" "GOTO" NEXT;
"IF" BB[K] > BBMAX "THEN" BBMAX:= BB[K]
"END";
"IF" K = M1 "THEN" M:= M1 "ELSE"
"BEGIN"
T:= D[M] - D[M1]; R:= BB[M1];
"IF" ABS(T) < TOL "THEN" S:= SQRT(R) "ELSE"
"BEGIN" W:= 2 / T; S:= W * R / (SQRT(W * W * R + 1) + 1)
"END"; "IF" K = M - 2 "THEN"
"BEGIN" D[M]:= D[M] + S; D[M1]:= D[M1] - S;
T:= - S / B[M1]; R:= SQRT(T * T + 1); COS:= 1 / R;
SIN:= T / R; ROTCOL(1,N,M1,M,A,COS,SIN); M:= M - 2
"END"
"ELSE"
"BEGIN" COUNT:= COUNT + 1;
"IF" COUNT > MAX "THEN" "GOTO" END;
LAMBDA:= D[M] + S; "IF" ABS(T) < TOL "THEN"
"BEGIN" W:= D[M1] - S;
"IF" ABS(W) < ABS(LAMBDA) "THEN" LAMBDA:= W
"END";
K:= K + 1; T:= D[K] - LAMBDA; COS:= 1; W:= B[K];
P:= SQRT(T * T + W * W); J1:= K;
"FOR" J:= K + 1 "STEP" 1 "UNTIL" M "DO"
"BEGIN" OLDCOS:= COS; COS:= T / P; SIN:= W / P;
DK1:= D[J] - LAMBDA; T:= OLDCOS * T;
D[J1]:= (T + DK1) * SIN * SIN + LAMBDA + T;
T:= COS * DK1 - SIN * W * OLDCOS; W:= B[J];
P:= SQRT(T * T + W * W); G:= B[J1]:= SIN * P;
BB[J1]:= G * G; ROTCOL(1, N, J1, J, A, COS, SIN);
J1:= J
"END";
D[M]:= COS * T + LAMBDA; "IF" T < 0 "THEN" B[M1]:= - G
"END" QRSTEP
"END";
"IF" M > 0 "THEN" "GOTO" IN;
END: EM[3]:= SQRT(BBMAX); EM[5]:= COUNT; QRISYMTRI:= M
"END" QRISYMTRI;
"EOP"
1SECTION 3.3.1.1.2 (JULY 1974) PAGE 1
AUTHORS: T.J.DEKKER AND W.HOFFMANN.
CONTRIBUTORS: W.HOFFMANN, J.G.VERWER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 730924.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS SEVEN PROCEDURES.
A) EIGVALSYM1 AND EIGVALSYM2 CALCULATE ALL EIGENVALUES, OR SOME
CONSECUTIVE EIGENVALUES INCLUDING THE LARGEST, OF A SYMMETRIC
MATRIX USING LINEAR INTERPOLATION ON A FUNCTION DERIVED FROM A
STURM SEQUENCE,
B) EIGSYM1 AND EIGSYM2 CALCULATE THE CORRESPONDING EIGENVECTORS
AS WELL, BY MEANS OF INVERSE ITERATION,
C) QRIVALSYM1 AND QRIVALSYM2 CALCULATE ALL EIGENVALUES OF A
SYMMETRIC MATRIX BY MEANS OF QR ITERATION,
D) QRISYM CALCULATES ALL EIGENVECTORS AS WELL IN THE SAME ITERATION
PROCESS.
EIGVALSYM1, EIGSYM1 AND QRIVALSYM1 USE IONAL ARRAY FOR
THE GIVEN SYMMETRIC MATRIX; THE OTHER PROCEDURES EXPECT THE MATRIX
TO BE STORED IN "ARRAY".
QRISYM DELIVERS THE EIGENVECTORS IN THE ARRAY THAT WAS USED FOR THE
ORIGINAL MATRIX IN CONTRAST WITH EIGSYM1 AND EIGSYM2 WHICH DELIVER
THE EIGENVECTORS IN AN EXTRA ARRAY.
WHEN ALL EIGENVALUES HAVE TO BE CALCULATED, THE PROCEDURES USING QR
ITERATION ARE PREFERABLE WITH RESPECT TO THEIR RUNNING TIME. WHEN
ALSO THE EIGENVECTORS HAVE TO BE CALCULATED THE PROCEDURES USING
INVERSE ITERATION ARE FASTER; HOWEVER, THE ONE USING QR ITERATION
USES LESS MEMORY SPACE.
KEYWORDS:
EIGENVALUES,
EIGENVECTORS,
SYMMETRIC MATRIX,
STURM-SEQUENCE,
INVERSE ITERATION,
QR ITERATION.
1SECTION 3.3.1.1.2 (JULY 1974) PAGE 2
SUBSECTION: EIGVALSYM2.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" EIGVALSYM2(A, N, NUMVAL, VAL, EM);
"VALUE" N, NUMVAL; "INTEGER" N, NUMVAL; "ARRAY" A, VAL, EM;
"CODE" 34153;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
NUMVAL: <ARITHMETIC EXPRESSION>;
THE SERIAL NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
A: <ARRAY IDENTIFIER>;
"ARRAY" A[1:N,1:N];
ENTRY: THE UPPER TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
GIVEN IN THE UPPER TRIANGULAR PART OF A (THE
ELEMENTS A[I,J], I<= J);
EXIT: THE DATA FOR HOUSEHOLDER'S BACK TRANSFORMATION
(WHICH ISN'T USED BY THIS PROCEDURE) IS DELIVERED
IN THE UPPER TRIANGULAR PART OF A;
THE ELEMENTS A[I,J] FOR I > J ARE NEITHER USED NOR CHANGED;
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[1:NUMVAL];
EXIT: THE NUMVAL LARGEST EIGENVALUES IN MONOTONICALLY
NON-INCREASING ORDER;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:3];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES;
EXIT: EM[1], THE INFINITY NORM OF THE ORIGINAL MATRIX,
EM[3], THE NUMBER OF ITERATIONS USED FOR
CALCULATING THE NUMVAL EIGENVALUES.
PROCEDURES USED:
TFMSYMTRI2 = CP34140,
VALSYMTRI = CP34151.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH:
THREE ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE USED.
RUNNING TIME:
ROUGHLY PROPORTIONAL TO N CUBED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE BODY OF EIGVALSYM2 CONSISTS OF TWO PROCEDURE STATEMENTS; THE
FIRST IS A CALL OF TFMSYMTRI2 TO TRANSFORM THE SYMMETRIC MATRIX
INTO A SIMILAR TRIDIAGONAL MATRIX BY MEANS OF HOUSEHOLDER'S
TRANSFORMATION; THE SECOND IS A CALL OF VALSYMTRI TO CALCULATE THE
DESIRED EIGENVALUES. OPERATION DETAILS OF BOTH PROCEDURES ARE GIVEN
IN THEIR DESCRIPTION.
1SECTION 3.3.1.1.2 (JULY 1974) PAGE 3
SUBSECTION: EIGVALSYM1.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" EIGVALSYM1(A, N, NUMVAL, VAL, EM);
"VALUE" N, NUMVAL; "INTEGER" N, NUMVAL; "ARRAY" A, VAL, EM;
"CODE" 34155;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
NUMVAL: <ARITHMETIC EXPRESSION>;
THE SERIAL NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
A: <ARRAY IDENTIFIER>;
"ARRAY" A[1:(N+1)*N//2];
ENTRY: THE UPPER TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
GIVEN IN SUCH A WAY THAT THE (I,J)-TH ELEMENT OF
THE MATRIX IS A[(J-1)*J//2+I], 1 <= I <= J <= N;
EXIT: THE DATA FOR HOUSEHOLDER'S BACK TRANSFORMATION
(WHICH ISN'T USED BY THIS PROCEDURE).
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[1:NUMVAL];
EXIT: THE NUMVAL LARGEST EIGENVALUES IN MONOTONICALLY
NON-INCREASING ORDER;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:3];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES;
EXIT: EM[1], THE INFINITY NORM OF THE ORIGINAL MATRIX,
EM[3], THE NUMBER OF ITERATIONS USED FOR
CALCULATING THE NUMVAL EIGENVALUES.
PROCEDURES USED:
TFMSYMTRI1 = CP34143,
VALSYMTRI = CP34151.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH:
THREE ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE USED.
RUNNING TIME:
ROUGHLY PROPORTIONAL TO N CUBED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE BODY OF EIGVALSYM1 CONSISTS OF TWO PROCEDURE STATEMENTS; THE
FIRST IS A CALL OF TFMSYMTRI1 TO TRANSFORM THE SYMMETRIC MATRIX
INTO A SIMILAR TRIDIAGONAL MATRIX BY MEANS OF HOUSEHOLDER'S
TRANSFORMATION; THE SECOND IS A CALL OF VALSYMTRI TO CALCULATE THE
DESIRED EIGENVALUES. OPERATION DETAILS OF BOTH PROCEDURES ARE GIVEN
IN THEIR DESCRIPTION.
1SECTION 3.3.1.1.2 (JULY 1974) PAGE 4
SUBSECTION: EIGSYM2.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" EIGSYM2(A, N, NUMVAL, VAL, VEC, EM);
"VALUE" N, NUMVAL; "INTEGER" N, NUMVAL; "ARRAY" A, VAL, VEC, EM;
"CODE" 34154;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
NUMVAL: <ARITHMETIC EXPRESSION>;
THE SERIAL NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
A: <ARRAY IDENTIFIER>;
"ARRAY" A[1:N,1:N];
ENTRY: THE UPPER TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
GIVEN IN THE UPPER TRIANGULAR PART OF A (THE
ELEMENTS A[I,J], I<= J);
EXIT: THE DATA FOR HOUSEHOLDER'S BACK TRANSFORMATION
IS DELIVERED IN THE UPPER TRIANGULAR PART OF A;
THE ELEMENTS A[I,J] FOR I > J ARE NEITHER USED NOR CHANGED;
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[1:NUMVAL];
EXIT: THE NUMVAL LARGEST EIGENVALUES IN MONOTONICALLY
NON-INCREASING ORDER;
VEC: <ARRAY IDENTIFIER>;
"ARRAY" VEC[1:N,1:NUMVAL];
EXIT: THE NUMVAL CALCULATED EIGENVECTORS, STORED COLUMN-
WISE, CORRESPONDING TO THE CALCULATED EIGENVALUES;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:9];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES,
EM[4], THE ORTHOGONALISATION PARAMETER (SEE METHOD
AND PERFORMANCE),
EM[6], THE TOLERANCE FOR THE EIGENVECTORS,
EM[8], THE MAXIMUM NUMBER OF INVERSE ITERATIONS
ALLOWED FOR THE CALCULATION OF EACH EIGEN-
VECTOR;
EXIT: EM[1], THE INFINITY NORM OF THE MATRIX,
EM[3], THE NUMBER OF ITERATIONS USED FOR
CALCULATING THE NUMVAL EIGENVALUES,
EM[5], THE NUMBER OF EIGENVECTORS INVOLVED IN THE
LAST GRAM-SCHMIDT ORTHOGONALISATION,
EM[7], THE MAXIMUM EUCLIDEAN NORM OF THE RESIDUES
OF THE CALCULATED EIGENVECTORS,
EM[9], THE LARGEST NUMBER OF INVERSE ITERATIONS
PERFORMED FOR THE CALCULATION OF SOME EIGEN-
VECTOR.
1SECTION 3.3.1.1.2 (JULY 1974) PAGE 5
PROCEDURES USED:
TFMSYMTRI2 = CP34140,
VALSYMTRI = CP34151,
VECSYMTRI = CP34152,
BAKSYMTRI2 = CP34141.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH:
THREE ONE-DIMENSIONAL REAL ARRAYS OF LENGTH N ARE DECLARED;
MOREOVER, VECSYMTRI USES FIVE ONE-DIMENSIONAL REAL ARRAYS
OF LENGTH N AND ONE BOOLEAN ARRAY OF LENGTH N.
RUNNING TIME:
ROUGHLY PROPORTIONAL TO N CUBED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE BODY OF EIGSYM2 CONSISTS OF FOUR PROCEDURE STATEMENTS;
THE FIRST IS A CALL OF TFMSYMTRI2 TO TRANSFORM THE SYMMETRIC
MATRIX INTO A SIMILAR TRIDIAGONAL MATRIX BY MEANS OF
HOUSEHOLDERS TRANSFORMATION,
THE SECOND IS A CALL OF VALSYMTRI TO CALCULATE THE DESIRED
EIGENVALUES,
THE THIRD IS A CALL OF VECSYMTRI TO CALCULATE THE CORRESPONDING
EIGENVECTORS AND
THE FOURTH IS A CALL OF BAKSYMTRI2 TO PERFORM THE BACK
TRANSFORMATION.
THE PARAMETERS EM[5], EM[7] AND EM[9] ARE GIVEN ITS VALUE IN THE
PROCEDURE VECSYMTRI. FOR A POSSIBLY SUBSEQUENT CALL OF VECSYMTRI
THE VALUE OF EM[5] IS NEEDED. WHEN CONSECUTIVE EIGENVALUES ARE TOO
CLOSE TOGETHER, THE CORRESPONDING EIGENVECTORS ARE NOT NECESSARILY
DELIVERED ORTHOGONAL BY INVERSE ITERATION (THE METHOD WHICH IS USED
IN VECSYMTRI). THEREFORE GRAM-SCHMIDT ORTHOGONALISATION IS APPLIED
ON THE EIGENVECTORS WHEN THE DISTANCE BETWEEN TWO CONSECUTIVE
EIGENVALUES IS SMALLER THAN EM[4].
FOR FURTHER DETAILS ONE IS REFERRED TO THE SPECIFIC PROCEDURE
DESCRIPTIONS.
1SECTION 3.3.1.1.2 (JULY 1974) PAGE 6
SUBSECTION: EIGSYM1.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"PROCEDURE" EIGSYM1(A, N, NUMVAL, VAL, VEC, EM);
"VALUE" N, NUMVAL; "INTEGER" N, NUMVAL; "ARRAY" A, VAL, VEC, EM;
"CODE" 34156;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE GIVEN MATRIX;
NUMVAL: <ARITHMETIC EXPRESSION>;
THE SERIAL NUMBER OF THE LAST EIGENVALUE TO BE CALCULATED;
A: <ARRAY IDENTIFIER>;
"ARRAY" A[1:(N+1)*N//2];
ENTRY: THE UPPER TRIANGLE OF THE SYMMETRIC MATRIX MUST BE
GIVEN IN SUCH A WAY THAT THE (I,J)-TH ELEMENT OF
THE MATRIX IS A[(J-1)*J//2+I], 1 <= I <= J <= N;
EXIT: THE DATA FOR HOUSEHOLDER'S BACK TRANSFORMATION;
VAL: <ARRAY IDENTIFIER>;
"ARRAY" VAL[1:NUMVAL];
EXIT: THE NUMVAL LARGEST EIGENVALUES IN MONOTONICALLY
NON-INCREASING ORDER;
VEC: <ARRAY IDENTIFIER>;
"ARRAY" VEC[1:N,1:NUMVAL];
EXIT: THE NUMVAL CALCULATED EIGENVECTORS, STORED COLUMN-
WISE, CORRESPONDING TO THE CALCULATED EIGENVALUES;
EM: <ARRAY IDENTIFIER>;
"ARRAY" EM[0:9];
ENTRY: EM[0], THE MACHINE PRECISION,
EM[2], THE RELATIVE TOLERANCE FOR THE EIGENVALUES,
EM[4], THE ORTHOGONALISATION PARAMETER (SEE METHOD
AND PERFORMANCE),
EM[6], THE TOLERANCE FOR THE EIGENVECTORS,
EM[8], THE MAXIMUM NUMBER OF INVERSE ITERATIONS
ALLOWED FOR THE CALCULATION OF EACH EIGEN-
VECTOR;
EXIT: EM[1], THE INFINITY NORM OF THE MATRIX,
EM[3], THE NUMBER OF ITERATIONS USED FOR
CALCULATING THE NUMVAL EIGENVALUES,
EM[5], THE NUMBER OF EIGENVECTORS INVOLVED IN THE
LAST GRAM-SCHMIDT ORTHOGONALISATION,