-
Notifications
You must be signed in to change notification settings - Fork 2
/
tominv.f
2425 lines (2396 loc) · 80.8 KB
/
tominv.f
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
CM
C->>> -------------------------------------------------> ems_tom_inv <<<
subroutine ems_tom_inv(
& rt_cod, rp_cn, rp_lvl, g_tt_da,
& n_r, mx_n_c, n_c, n_a_el, n_inv_sing, fill_fac,
& vr_in_r,
& mtx_r_v, mtx_r_ix, mtx_c_sa,
& mx_n_eta, mx_n_eta_el, n_eta, n_eta_el, n_lo_eta,
& eta_v, eta_ix, eta_sa,
& lo_eta_pv_in_r, up_eta_pv_in_r,
& pv_tl, wr_eta_tl, c_rlv_tl, unit_v_tl,
& l_pv_rsdu_er_tl,
& wk_a, bs_c_at, r_k_a, pv_r_n_or_mrt,
& ck_l_pv_wk_a, ds, is)
implicit none
include 'EMSV.INC'
include 'TOMINV.INC'
CM IF (emsol_tt .EQ. 1) THEN
C? include 'EMSTT.INC'
CM ENDIF
integer rt_cod, rp_cn, rp_lvl, g_tt_da
integer n_r, n_c, mx_n_c, n_a_el, n_inv_sing
double precision fill_fac
integer mx_n_eta, mx_n_eta_el, n_eta, n_eta_el, n_lo_eta
integer mtx_r_ix(0:n_a_el)
integer mtx_c_sa(0:n_c + 1)
integer eta_ix(0:mx_n_eta_el)
integer eta_sa(0:mx_n_eta+1)
integer lo_eta_pv_in_r(0:n_r)
integer up_eta_pv_in_r(0:n_r)
integer is(0:*)
double precision pv_tl, wr_eta_tl, c_rlv_tl, unit_v_tl
double precision l_pv_rsdu_er_tl
double precision mtx_r_v(0:n_a_el)
double precision wk_a(0:n_r)
double precision eta_v(0:mx_n_eta_el)
double precision ck_l_pv_wk_a(0:n_r)
double precision ds(0:*)
integer bs_c_at(1:n_r)
integer r_k_a(1:n_r)
integer pv_r_n_or_mrt(1:n_r)
integer vr_in_r(0:n_r)
double precision mn_el_z, c_mx
double precision lc_v, abs_lc_v, c_rlv_mn, neg_rcp_pv_v
integer n_bs_non_ze
integer n_lg, n_in_bump, n_ab_bump, n_bw_bump
integer r_n, vr_n, i, k, c_n, el_n, inner_el_n
integer n_in_un_pv_r, pv_r_cdd
integer n_in_un_pv_r_m1, se_1_fh_at_sa
integer n_redo_c
integer lc_mtx_c_sa, lc_mtx_c_fh
integer mrt_k, stp, mn_k, mrt_k_dl
integer se_1_sa, se_1_fh, se_1_en_n
integer se_2_sa, se_2_fh, se_2_en_n
integer se_3_sa, se_3_fh, se_3_en_n
integer f_eta_n, f_eta_el_n
integer f_bmp_lo_eta_n, f_bmp_lo_el_n
integer fm_bmp_lo_eta_n
integer l_lo_eta_n, l_lo_el_n_p1, d_l_lo_el_n_p1
integer l_up_eta_n, l_up_el_n_m1, d_l_up_el_n_m1
integer l_eta_n, l_eta_el_n_p1, d_l_eta_el_n_p1
integer bst_eta_el_n, l_up_el_n, l_eta_el_n, n_eta_el_qy
integer eta_n, n_triang_ps
integer n_mv_c
logical pv_ok, mv_c
integer sus_ftran_n_op
integer sus_ftran_mx_n_op
integer dl_el_n
CM IF (emsol_da .EQ. 1) THEN
C?c integer bmp_dim
C?c integer n_ty_1_ftran
C?c integer n_ty_2_ftran
C?c integer n_ftran_eta
C?c integer n_ftran_eta_skp_n
C?c integer n_ftran_eta_skp_en
C?c integer su_n_ftran_eta
C?c integer n_no_apply_ftran_eta
C?c double precision pct
CM ENDIF
CM IF (emsol_deb .EQ. 1) THEN
C? integer lc_rt_cod
C?c logical er_fd
CM ENDIF
CM IF (emsol_tt .EQ. 1) THEN
C? logical g_tom_tt_lvl2
C? logical g_tom_tt_lvl3
CM ENDIF
c
c pv_r_n_or_mrt is used to store the pivot row for the columns in
c SET 2 (below the bump columns) and SET 3 (the logicals)
c and the merit for the columns in the bump (part of SET 1).
c The SET 3 storage is only used for a short time as temporary
C storage while the logicals are repositioned in vr_in_r.
c if (rp_cn .ge. 0) write(rp_cn, *) 'on entry vr_in_r='
c if (rp_cn .ge. 0) write(rp_cn, *) (vr_in_r(i), i=1, n_r)
c
c Tolerances [default---see ems_iz_rl_ct_vr_df]
c
c pv_tl [1d-6]: The pivot tolerance: Values less than this are
c written to the eta file but are not considered
c for pivots.
c
c wr_eta_tl [1d-9]: The drop tolerance: Values less than this are
c not written to the eta file
c
c c_rlv_tl [1d-2]: 'u' in threshold pivoting
c
c unit_v_tl [1d-6]: Used to test whether a value can be taken as 1.
c
c=======================================================================
c
c The value of tom_inv_sing_msk determines the tests and subsequent
c action in rtelation to the final pivot.
c
c if iand(tom_inv_sing_msk, *) then
c
c tom_inv_sing_tl_ck: a test is done if the pivot is small
c
c tom_inv_sing_al_ck: a test is always done
c
c tom_inv_sing_tl_alw: small pivots passing the test are allowed
c
c tom_inv_sing_al_alw: small pivots are always allowed
c
c NB
c tom_inv_sing_al_ck => tom_inv_sing_tl_ck
c tom_inv_sing_al_alw => tom_inv_sing_tl_alw
c
c Suggested values
c
c tom_inv_sing_msk = 0: Reverts to old Tomlin INVERT
c tom_inv_sing_msk = 1: Tests small pivots
c tom_inv_sing_msk = 1+2 = 3: Tests all pivots
c tom_inv_sing_msk = 1+4 = 5: Tests small and allows OK pivots
c tom_inv_sing_msk = 3+4 = 7: Tests all pivots and allows OK pivots
c tom_inv_sing_msk = 4+8 = 12: Allows all pivots without testing
c tom_inv_sing_msk = 1+12= 13: Allows all pivots, testing small
c tom_inv_sing_msk = 3+12= 15: Allows all pivots, testing all
c
c=======================================================================
CM IF (emsol_tt .EQ. 1) THEN
C? g_tom_tt_lvl2 = iand(g_tt_da, bt11) .ne. 0
C? g_tom_tt_lvl3 = iand(g_tt_da, bt12) .ne. 0
CM ENDIF
rt_cod = 0
n_eta = 0
n_eta_el = 0
n_ab_bump = 0
n_in_bump = 0
n_bw_bump = 0
n_inv_sing = 0
n_redo_c = 0
fill_fac = one
l_lo_eta_n = 0
f_eta_el_n = 1
l_eta_el_n_p1 = f_eta_el_n
f_eta_n = 1
l_eta_n = 0
se_1_sa = 1
se_1_fh = se_1_sa - 1
se_3_sa = n_r + 1
se_3_fh = se_3_sa - 1
c
c se_1_sa remains 1
c se_3_fh remains n_r
c
sus_ftran_mx_n_op = lo_eta_pv_in_r(0)
CM IF (emsol_da .EQ. 1) THEN
C?c n_ty_1_ftran = 0
C?c n_ty_2_ftran = 0
C?c n_ftran_eta_skp_n = 0
C?c n_ftran_eta_skp_en = 0
C?c su_n_ftran_eta = 0
CM ENDIF
c
c Pull out logicals for SET 3; rest in SET 1
c
do 110, r_n = 1, n_r
c
c Initialise lo_eta_pv_in_r (up_eta_pv_in_r) which contains the eta
c number of any L-eta (U-eta) which has a pivot in that row. Those
c rows in which no L-eta (U-eta) has a pivot must have a large eta
c number since the minimum eta number will be required to start
c FTRAN.
c
lo_eta_pv_in_r(r_n) = mx_n_eta+1
up_eta_pv_in_r(r_n) = mx_n_eta+1
if (vr_in_r(r_n).gt.mx_n_c) then
se_3_sa = se_3_sa - 1
pv_r_n_or_mrt(se_3_sa) = vr_in_r(r_n)
bs_c_at(se_3_sa) = vr_in_r(r_n)
else
se_1_fh = se_1_fh + 1
bs_c_at(se_1_fh) = vr_in_r(r_n)
end if
r_k_a(r_n) = -1
C vr_in_r(r_n) = 0
110 continue
c
c Set 2 is initially empty: It will be formed between sets 1 and 3!
c
se_2_fh = se_3_sa - 1
se_2_sa = se_3_sa
c
c Loop through the logicals:
c * setting vr_in_r so they are solved for in their home row;
c * setting r_k_a = 0 for each row with a basic logical.
c
do 210, se_3_en_n = se_3_sa, se_3_fh
vr_n = pv_r_n_or_mrt(se_3_en_n)
r_k_a(vr_n-mx_n_c) = 0
vr_in_r(vr_n-mx_n_c) = vr_n
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & ' Logical: ', vr_n-mx_n_c,
C? & ' Pivoting in row ', vr_n-mx_n_c
C? endif
CM ENDIF
210 continue
c
c The initialisation of bs_c_at and r_k_a is complete:
c
c bs_c_at(se_en_n) is structural for se_en_n = 1, se_1_fh
c bs_c_at(se_en_n) is logical for se_en_n = se_3_sa, se_3_fh
c
c r_k_a(r_n) = 0 if bs_c_at(r_n) is logical.
c r_k_a(r_n) = -1 if bs_c_at(r_n) is structural
c
c pv_r_n_or_mrt is partially initialised:
c
c pv_r_n_or_mrt(se_en_n) for se_en_n = se_3_sa, se_3_fh
c . is the logical variable number (NOT its
c . pivotal row number---why?)
c
c In 1 pass pull out some vectors below bump and get row counts
c
n_lg = se_3_fh - se_3_sa + 1
n_bs_non_ze = n_lg
n_triang_ps = 1
if (se_1_fh .eq. 0) go to 3700
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl2) call ems_tt_rec(tom_triang_tt, -1)
CM ENDIF
c
c Cycle through the set 1 entries.
c
se_1_en_n = se_1_sa
c
c ==================================================================
c Top of first triangularisation section
c
c For each column in set 1
c * Accumulate the number of nonzeros in the basis---not used
c * Count the number of entries in rows for which a pivot has not
c . yet been decided (rows with basic logical) and update the row
c . count for rows with basic structural.
c . o If there is exactly one entry in non-pivotal rows then the
c . column can be moved into set 2 with a pivot in this row.
c . o If there is more than one entry in non-pivotal rows then the
c . column remains in set 1.
c . o If there is no entry in non-pivotal rows then the matrix is
c . singular and the column is flagged for removal.
c
300 continue
c_n = bs_c_at(se_1_en_n)
n_in_un_pv_r = 0
n_bs_non_ze = n_bs_non_ze + (mtx_c_sa(c_n+1)-mtx_c_sa(c_n))
do 310, el_n = mtx_c_sa(c_n), mtx_c_sa(c_n+1)-1
r_n = mtx_r_ix(el_n)
if (r_k_a(r_n) .ge. 0) go to 310
n_in_un_pv_r = n_in_un_pv_r + 1
r_k_a(r_n) = r_k_a(r_n) - 1
pv_r_cdd = r_n
310 continue
n_in_un_pv_r_m1 = n_in_un_pv_r - 1
if (n_in_un_pv_r_m1 .eq. 0) then
c
c Column singleton:
c Exactly one nonzero in rows which have not yet been pivoted on.
c Move entry into set 2 by swapping it with the last entry in set 1
c and shifting the finish of set 1 and start of set 2 by 1.
c
bs_c_at(se_1_en_n) = bs_c_at(se_1_fh)
se_1_fh = se_1_fh - 1
se_2_sa = se_2_sa - 1
bs_c_at(se_2_sa) = c_n
c
c Indicate that this set entry has been pivoted on by recording its
c pivotal row number, zeroing its row count and that it is solved
c for in this row.
c
pv_r_n_or_mrt(se_2_sa) = pv_r_cdd
r_k_a(pv_r_cdd) = 0
vr_in_r(pv_r_cdd) = c_n
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & ' Column singleton: ', c_n,
C? & ' Pivoting in row ', pv_r_cdd
C? endif
CM ENDIF
c
c Go back to consider next column in set 1---the one which was
c swapped in to the current set 1 entry---unless set 1 is finished
c (in which case the column had been at the end of set 1 entry and
c was swapped with itself).
c
if (se_1_en_n .gt. se_1_fh) go to 400
go to 300
else if (n_in_un_pv_r_m1 .gt. 0) then
c
c More than one nonzero in rows which have not yet been pivoted on.
c Column remains in set 1. Go back to consider next column in set 1
c unless set 1 is finished.
c
if (se_1_en_n .ge. se_1_fh) go to 400
se_1_en_n = se_1_en_n + 1
go to 300
else
c
c No nonzero in rows which have not yet been pivoted on.
c Matrix is singular so mark the column as one to be removed.
c Column is removed from set 1 by swapping it with the last entry in
c set 1 and shifting the finish of set 1.
c
n_inv_sing = n_inv_sing + 1
call ems_se_iv_ty(1, rp_cn, c_n, ds, is)
bs_c_at(se_1_en_n) = bs_c_at(se_1_fh)
se_1_fh = se_1_fh - 1
if (se_1_en_n .gt. se_1_fh) go to 400
go to 300
end if
400 continue
c
c Set 2 now contains structurals which can be pivoted on after the
c bump and before the logicals.
c
n_bw_bump = se_2_fh - se_2_sa + 1
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & 'Triang', n_triang_ps, ' n_ab_bump =', n_ab_bump,
C? & ' n_bw_bump =', n_bw_bump
C? endif
CM ENDIF
c
c Bottom of first triangularisation section
c ==================================================================
c Top of second and subsequent triangularisation section
c
c Before starting the second and subsequent triangularisation pass,
c see whether INVERT has finished.
c
410 continue
if (se_1_fh .eq. 0) then
c
c All structurals are below the bump (ie there is no bump!) This
c happens in network problems.
c
n_bw_bump = se_2_fh - se_2_sa + 1
c
c The triangularisation phase is finished.
c
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl2) call ems_tt_rec(-tom_triang_tt, -1)
CM ENDIF
eta_sa(l_eta_n+1) = l_eta_el_n_p1
c
c Record the (current) last lower triangular eta number and element
c number+1
c
l_lo_eta_n = l_eta_n
l_lo_el_n_p1 = l_eta_el_n_p1
c
c Write in the etas for columns below the bump.
c
go to 3700
end if
c
c Have to do the triangularisation pass so increment the number of
c passes.
c
n_triang_ps = n_triang_ps + 1
c
c Record the finish of set 1 at the start of this pass and start
c with the first entry in set 1.
c
se_1_fh_at_sa = se_1_fh
se_1_en_n = se_1_sa
420 continue
c_n = bs_c_at(se_1_en_n)
n_in_un_pv_r = 0
lc_mtx_c_sa = mtx_c_sa(c_n)
lc_mtx_c_fh = mtx_c_sa(c_n+1)-1
el_n = lc_mtx_c_sa
go to 510
c
c Look for an entry in the column which is in a singleton row
c ---a row for which r_k_a(r_n) = -2 (since r_k_a(r_n) is
c initialised to -1)
c
500 continue
if (el_n .ge. lc_mtx_c_fh) go to 600
el_n = el_n + 1
510 continue
r_n = mtx_r_ix(el_n)
if (r_k_a(r_n) .ge. 0) go to 500
if (r_k_a(r_n) .eq. -2) then
c
c Row singleton:
c The column will appear before the bump so its L-eta can be
c written now.
c
n_ab_bump = n_ab_bump + 1
c
c Check that there is enough space to write the eta
c
i = lc_mtx_c_fh - lc_mtx_c_sa
if (l_eta_el_n_p1 + i .gt. mx_n_eta_el) go to 4000
c
c Increment the number of the last eta (in the file to date),
c store the start of this eta, the pivotal row and the negative
c reciprocal of the pivot, and then increment the record of the
c number of the last eta element (in the file to date).
c
l_eta_n = l_eta_n + 1
eta_sa(l_eta_n) = l_eta_el_n_p1
eta_ix(l_eta_el_n_p1) = r_n
eta_v(l_eta_el_n_p1) = -1/mtx_r_v(el_n)
l_eta_el_n_p1 = l_eta_el_n_p1 + 1
lo_eta_pv_in_r(r_n) = l_eta_n
c
c Set the row count to 1---why?
c
r_k_a(r_n) = 1
vr_in_r(r_n) = c_n
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & ' Row singleton: ', c_n,
C? & ' Pivoting in row ', r_n
C? endif
CM ENDIF
el_n = el_n - 1
520 continue
c
c Store the entries in the eta either
c * up to the pivot---in the first pass
c * or after the pivot---in the second pass when
c . lc_mtx_c_sa = el_n + 2 and el_n = lc_mtx_c_fh.
c
do 530, inner_el_n = lc_mtx_c_sa, el_n
r_n = mtx_r_ix(inner_el_n)
c
c If the row is non-pivotal, update the count of entries in
c columns which may form the bump since the current column
c is now before the bump.
c
if (r_k_a(r_n) .lt. 0) r_k_a(r_n) = r_k_a(r_n) + 1
eta_ix(l_eta_el_n_p1) = r_n
eta_v(l_eta_el_n_p1) = mtx_r_v(inner_el_n)
l_eta_el_n_p1 = l_eta_el_n_p1 + 1
530 continue
if (el_n .lt. lc_mtx_c_fh) then
c
c If there are still entries to store then change the loop limits
c and make a second pass.
c
lc_mtx_c_sa = el_n + 2
el_n = lc_mtx_c_fh
go to 520
end if
c
c Remove the column from set 1 by swapping it with the last entry
c in set 1 and decreasing the finish of set 1.
c Jump to a line which is just a jump back to start with the next
c set 1 entry---if there are any more.
c
bs_c_at(se_1_en_n) = bs_c_at(se_1_fh)
se_1_fh = se_1_fh - 1
go to 700
end if
c
c Increment the count of entries in this column in non-pivotal rows
c which, if equal to 0 (1) indicates singularity (column singleton).
c
n_in_un_pv_r = n_in_un_pv_r + 1
pv_r_cdd = r_n
c
c Jump back to consider the next entry in the column.
c
go to 500
c
c All the entries in the column have been considered: none is in
c a singleton row so see whether this column is a singleton.
c
600 continue
n_in_un_pv_r_m1 = n_in_un_pv_r - 1
if (n_in_un_pv_r_m1 .eq. 0) then
c
c Column singleton:
c Exactly one nonzero in rows which have not yet been pivoted on.
c Move entry into set 2 by swapping it with the last entry in set 1
c and shifting the finish of set 1 and start of set 2 by 1.
c
bs_c_at(se_1_en_n) = bs_c_at(se_1_fh)
se_1_fh = se_1_fh - 1
se_2_sa = se_2_sa - 1
bs_c_at(se_2_sa) = c_n
c
c Indicate that this set entry has been pivoted on by recording its
c pivotal row number, zeroing its row count and that it is solved
c for in this row.
c
pv_r_n_or_mrt(se_2_sa) = pv_r_cdd
r_k_a(pv_r_cdd) = 0
vr_in_r(pv_r_cdd) = c_n
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & ' Column singleton: ', c_n,
C? & ' Pivoting in row ', pv_r_cdd
C? endif
CM ENDIF
go to 700
else if (n_in_un_pv_r_m1 .gt. 0) then
c
c More than one nonzero in rows which have not yet been pivoted on.
c Column remains in set 1. Go back to consider next column in set 1
c unless set 1 is finished.
c
if (se_1_en_n .ge. se_1_fh) go to 710
se_1_en_n = se_1_en_n + 1
go to 420
else
c
c No nonzero in rows which have not yet been pivoted on.
c Matrix is singular so mark the column as one to be removed.
c Column is removed from set 1 by swapping it with the last entry in
c set 1 and shifting the finish of set 1.
c
n_inv_sing = n_inv_sing + 1
call ems_se_iv_ty(2, rp_cn, c_n, ds, is)
bs_c_at(se_1_en_n) = bs_c_at(se_1_fh)
se_1_fh = se_1_fh - 1
end if
700 continue
c
c If not all the columns in set 1 have been considered then jump
c back to consider the next
c
if (se_1_en_n .le. se_1_fh) go to 420
c
c Jump to here if set 1 is known to have been finished when the
c column is known to be a non-singleton.
c Saves just an increase in se_1_en_n to se_1_fh+1.
c
710 continue
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & 'Triang', n_triang_ps, ' n_ab_bump =', n_ab_bump,
C? & ' n_bw_bump =', se_2_fh - se_2_sa + 1
C? endif
CM ENDIF
c
c If the number of columns in set 1 has been reduced on this pass,
c repeat the triangularisation pass.
c
if (se_1_fh .ne. se_1_fh_at_sa) go to 410
c
c Store the end of the last eta (in the file to date).
c
eta_sa(l_eta_n+1) = l_eta_el_n_p1
c
c The triangularisation phase is finished.
c
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl2) call ems_tt_rec(-tom_triang_tt, -1)
CM ENDIF
c
c Bottom of second and subsequent triangularisation section.
c ==================================================================
c The bump is now finalised:
c Set 1 contains the structurals in the bump
c Set 2 contains the structurals after the bump (column singletons).
c Set 3 contains the logicals
c
c Between sets 1 and 2 are the structurals which are before the bump
c (row singletons) and singular structurals.
c
n_in_bump = se_1_fh - se_1_sa + 1
n_bw_bump = se_2_fh - se_2_sa + 1
c
c Record the (current) last lower triangular eta number and element
c number+1
c
l_lo_eta_n = l_eta_n
l_lo_el_n_p1 = l_eta_el_n_p1
c
c There is no bump so write in the etas for columns below the bump.
c
if (se_1_fh .eq. 0) go to 3700
c
c Get column merit counts for set 1: this is the sum of counts of
c the non-pivotal rows in each column.
c
do 820, se_1_en_n = se_1_sa, se_1_fh
c_n = bs_c_at(se_1_en_n)
mrt_k = 0
do 810, el_n = mtx_c_sa(c_n), mtx_c_sa(c_n+1)-1
r_n = mtx_r_ix(el_n)
if (r_k_a(r_n) .lt. 0) mrt_k = mrt_k - (r_k_a(r_n)+1)
810 continue
pv_r_n_or_mrt(se_1_en_n) = mrt_k
820 continue
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 2) then
C? do 830, se_1_en_n = se_1_sa, se_1_fh
C? if (rp_cn .ge. 0) write(rp_cn, 9110)'before_sort = ',
C? & bs_c_at(se_1_en_n), pv_r_n_or_mrt(se_1_en_n)
C? 830 continue
C? endif
CM ENDIF
c =================================================================
c Top of shell sort of set 1 columns by increasing merit count.
c
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl2) call ems_tt_rec(tom_srt_tt, -1)
CM ENDIF
stp = 1
910 continue
stp = 3*stp + 1
if (se_1_fh .ge. stp) go to 910
do 940, k = 1, i_inf
stp = stp/3
if (stp .lt. 1) go to 950
do 930, i = stp+1, se_1_fh
mrt_k = pv_r_n_or_mrt(i)
if (pv_r_n_or_mrt(i-stp) .le. mrt_k) go to 930
c_n = bs_c_at(i)
se_1_en_n = i
920 continue
pv_r_n_or_mrt(se_1_en_n) = pv_r_n_or_mrt(se_1_en_n-stp)
bs_c_at(se_1_en_n) = bs_c_at(se_1_en_n-stp)
se_1_en_n = se_1_en_n - stp
if (se_1_en_n .gt. stp) then
if (pv_r_n_or_mrt(se_1_en_n-stp) .gt. mrt_k) go to 920
end if
pv_r_n_or_mrt(se_1_en_n) = mrt_k
bs_c_at(se_1_en_n) = c_n
930 continue
940 continue
950 continue
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 2) then
C? do 960, se_1_en_n = se_1_sa, se_1_fh
C? if (rp_cn .ge. 0) write(rp_cn, 9110)'after_sort =',
C? & bs_c_at(se_1_en_n), pv_r_n_or_mrt(se_1_en_n)
C? 960 continue
C? endif
CM ENDIF
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl2) call ems_tt_rec(-tom_srt_tt, -1)
CM ENDIF
c
c Bottom of shell sort of set 1 columns by increasing merit count.
c =================================================================
c
c Set up for writing bump etas; below bump not yet written out.
c Bump lo elements written to continuation of lo part of file.
c Bump up elements written to end of file. Temporary values assigned
c to l_eta_n, to use FTRAN as half-tran. Proper
c values later copied back.
c
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl2) call ems_tt_rec(tom_bmp_inv_tt, -1)
CM ENDIF
c f_bmp_lo_eta_n = l_eta_n + 1
c f_bmp_lo_el_n = l_eta_el_n_p1
f_bmp_lo_eta_n = l_lo_eta_n + 1
f_bmp_lo_el_n = l_lo_el_n_p1
c
c Zero the work vector
c
do 1010, r_n = 1, n_r
wk_a(r_n) = 0
1010 continue
c
c n_eta_el_qy is ?
c
n_eta_el_qy = 0
l_up_el_n_m1 = mx_n_eta_el
l_up_eta_n = mx_n_eta + 1
eta_sa(mx_n_eta+1) = mx_n_eta_el+1
c
c Loop over the sorted columns in the bump.
c
n_mv_c = 0
CM IF (emsol_da .EQ. 1) THEN
C?c bmp_dim = se_1_fh - se_1_sa + 1
CM ENDIF
do 1500, se_1_en_n = se_1_sa, se_1_fh
1100 continue
c_n = bs_c_at(se_1_en_n)
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9000)
C? & '1: Set 2 entry ', se_1_en_n, ' Column ', c_n
C? endif
CM ENDIF
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9110)'after sort =', c_n
C? endif
CM ENDIF
d_l_up_el_n_m1 = l_up_el_n_m1
d_l_lo_el_n_p1 = l_lo_el_n_p1
c
c Check that there is space for n_r+1 eta elements.
c
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & ' d_l_up_el_n_m1 =', d_l_up_el_n_m1,
C? & ' d_l_lo_el_n_p1 =', d_l_lo_el_n_p1,
C? & ' n_r =', n_r
C? endif
CM ENDIF
if (d_l_up_el_n_m1 - d_l_lo_el_n_p1 .lt. n_r) then
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl2) call ems_tt_rec(tom_bmp_inv_tt, -1)
CM ENDIF
go to 4000
endif
c
c ==================================================================
c Top of FIRST unpack matrix column section
c
c Unpack the matrix column, storing its values in the work array and
c each row index in the U-eta or L-eta, according to whether the row
c is pivotal or not.
c
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl3) call ems_tt_rec(tom_unpk_mtx_c1_tt, -1)
CM ENDIF
fm_bmp_lo_eta_n = l_lo_eta_n+1
c i = f_bmp_lo_el_n
do 1110, el_n = mtx_c_sa(c_n), mtx_c_sa(c_n+1)-1
r_n = mtx_r_ix(el_n)
wk_a(r_n) = mtx_r_v(el_n)
fm_bmp_lo_eta_n = min(lo_eta_pv_in_r(r_n), fm_bmp_lo_eta_n)
if (r_k_a(r_n) .ge. 0) then
c
c Row is pivotal so this entry is part of the U-eta.
c
if (d_l_up_el_n_m1 .lt. 0) then
CM IF (emsol_deb .EQ. 1) THEN
C? call ems_baso(lc_rt_cod, ds, 26, 1)
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9899)
C? & eta_ix(d_l_up_el_n_m1)
C? endif
CM ENDIF
rt_cod = ior(rt_cod, tom_inv_er_bt)
endif
eta_ix(d_l_up_el_n_m1) = r_n
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & 'unpk_upper i, eta_ix(i), wk_a(eta_ix(i))',
C? & d_l_up_el_n_m1, r_n, wk_a(r_n)
C? endif
CM ENDIF
d_l_up_el_n_m1 = d_l_up_el_n_m1 - 1
go to 1110
end if
c
c Row is non-pivotal so this is entry is part of the L-eta or pivot
c
eta_ix(d_l_lo_el_n_p1) = r_n
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & 'unpk_lower i, r_n, wk_a(r_n)',
C? & d_l_lo_el_n_p1, r_n, wk_a(r_n)
C? endif
CM ENDIF
d_l_lo_el_n_p1 = d_l_lo_el_n_p1 + 1
1110 continue
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl3) call ems_tt_rec(-tom_unpk_mtx_c1_tt, -1)
CM ENDIF
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 0) then
C? if (rp_cn .ge. 0) write(rp_cn, 9100)
C? & ' d_l_up_el_n_m1 =', d_l_up_el_n_m1,
C? & ' d_l_lo_el_n_p1 =', d_l_lo_el_n_p1
C? endif
CM ENDIF
c
c End of FIRST unpack matrix column section
c ==================================================================
c
c ==================================================================
c Start of FIRST FTRAN section
c
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl3) call ems_tt_rec(tom_ftran1_tt, -1)
CM ENDIF
n_eta_el_qy = 4*n_eta_el_qy
CM IF (emsol_da .EQ. 1) THEN
C?c n_ty_1_ftran = n_ty_1_ftran + 1
CM ENDIF
if (sus_ftran_mx_n_op .lt. 0) then
c
c Don't do any super-sparse FTRAN operations
c
fm_bmp_lo_eta_n = f_bmp_lo_eta_n
goto 1250
else
c
c Jump to the end of FTRAN if no operations need be applied.
c
if (fm_bmp_lo_eta_n .gt. l_lo_eta_n) goto 1275
c
c Skip the super-sparse FTRAN operation if none is to be applied
c (ie all that is being done is modify the first parameter of the
c normal loop).
c
if (sus_ftran_mx_n_op .eq. 0) goto 1250
endif
c
c Apply the first sus_ftran_mx_n_op etas by exploiting
c lo_eta_pv_in_r to skip etas
c
sus_ftran_n_op = 0
1200 continue
sus_ftran_n_op = sus_ftran_n_op + 1
CM IF (emsol_da .EQ. 1) THEN
C?c
C?c Decrementing n_ftran_eta_skp_n is necessary since all etas before
C?c the ultimate value of fm_bmp_lo_eta_n are assumed to be skipped
C?c in the accounting below 1275
C?c
C?c n_ftran_eta_skp_n = n_ftran_eta_skp_n - 1
CM ENDIF
eta_n = fm_bmp_lo_eta_n
el_n = eta_sa(eta_n)
r_n = eta_ix(el_n)
c
c Better to hit the occasional zero in the RHS than have to test for
c zero every time in loop 1220
c
if (wk_a(r_n) .eq. zero) goto 1210
lc_v = wk_a(r_n)*eta_v(el_n)
wk_a(r_n) = lc_v
n_eta_el_qy =
& n_eta_el_qy + eta_sa(eta_n+1)-1 - eta_sa(eta_n)
do 1205, el_n = el_n+1, eta_sa(eta_n+1)-1
r_n = eta_ix(el_n)
if (wk_a(r_n) .eq. zero) then
wk_a(r_n) = lc_v*eta_v(el_n)
if (r_k_a(r_n) .ge. 0) then
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 1) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & ' New Upper ',
C? & d_l_up_el_n_m1, r_n, wk_a(r_n)
C? endif
CM ENDIF
eta_ix(d_l_up_el_n_m1) = r_n
d_l_up_el_n_m1 = d_l_up_el_n_m1 - 1
go to 1205
end if
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 1) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & ' New Lower ',
C? & d_l_lo_el_n_p1, r_n, wk_a(r_n)
C? endif
CM ENDIF
eta_ix(d_l_lo_el_n_p1) = r_n
d_l_lo_el_n_p1 = d_l_lo_el_n_p1 + 1
go to 1205
end if
wk_a(r_n) = wk_a(r_n) + lc_v*eta_v(el_n)
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 2) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & ' Updated value ',
C? & d_l_lo_el_n_p1, r_n, wk_a(r_n)
C? endif
CM ENDIF
if (wk_a(r_n) .ne. zero) then
if (lo_eta_pv_in_r(r_n) .gt. eta_n) fm_bmp_lo_eta_n =
& min(lo_eta_pv_in_r(r_n), fm_bmp_lo_eta_n)
endif
1205 continue
1210 continue
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl3) call ems_tt_rec(tom_ftran_eta_srch_tt, -1)
CM ENDIF
fm_bmp_lo_eta_n = l_lo_eta_n+1
c
c Go through the entries which are (potentially) in the U-eta,
c looking to see which eta (if any) needs to be applied next
c No entries which are (potentially) in the L-eta can have a pivot
c in that row.
c
do 1220, el_n = l_up_el_n_m1, d_l_up_el_n_m1+1, -1
r_n = eta_ix(el_n)
if (lo_eta_pv_in_r(r_n) .gt. eta_n) fm_bmp_lo_eta_n =
& min(lo_eta_pv_in_r(r_n), fm_bmp_lo_eta_n)
1220 continue
CM IF (emsol_tt .EQ. 1) THEN
C? if (g_tom_tt_lvl3) call ems_tt_rec(-tom_ftran_eta_srch_tt, -1)
CM ENDIF
c
c Jump to the end of FTRAN if no further operations need be applied.
c
if (fm_bmp_lo_eta_n .gt. l_lo_eta_n) goto 1275
c
c Jump to the normal loop if the limit on the number of super-sparse
c FTRAN operations has been reached
c
if (sus_ftran_n_op .lt. sus_ftran_mx_n_op) goto 1200
1250 continue
do 1260, eta_n = fm_bmp_lo_eta_n, l_lo_eta_n
el_n = eta_sa(eta_n)
r_n = eta_ix(el_n)
CM IF (emsol_da .EQ. 1) THEN
C? if (wk_a(r_n) .eq. zero) then
C?c n_ftran_eta_skp_en = n_ftran_eta_skp_en + 1
C? goto 1260
C? endif
CM ELSE
if (wk_a(r_n) .eq. zero) go to 1260
CM ENDIF
lc_v = wk_a(r_n)*eta_v(el_n)
wk_a(r_n) = lc_v
n_eta_el_qy =
& n_eta_el_qy + eta_sa(eta_n+1)-1 - eta_sa(eta_n)
do 1255, el_n = el_n + 1, eta_sa(eta_n+1)-1
r_n = eta_ix(el_n)
if (wk_a(r_n) .eq. zero) then
c
c new element
c
wk_a(r_n) = lc_v*eta_v(el_n)
if (r_k_a(r_n) .ge. 0) then
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 1) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & ' New Upper ',
C? & d_l_up_el_n_m1, r_n, wk_a(r_n)
C? endif
CM ENDIF
c
c part of upper
c
eta_ix(d_l_up_el_n_m1) = r_n
d_l_up_el_n_m1 = d_l_up_el_n_m1 - 1
go to 1255
end if
c
c pivot element or part of lower
c
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 1) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & ' New Lower ',
C? & d_l_lo_el_n_p1, r_n, wk_a(r_n)
C? endif
CM ENDIF
eta_ix(d_l_lo_el_n_p1) = r_n
d_l_lo_el_n_p1 = d_l_lo_el_n_p1 + 1
go to 1255
end if
c
c element previously recorded
c
wk_a(r_n) = wk_a(r_n) + lc_v*eta_v(el_n)
CM IF (emsol_deb .EQ. 1) THEN
C? if (rp_lvl .gt. 2) then
C? if (rp_cn .ge. 0) write(rp_cn, 9220)
C? & ' Updated value ',
C? & d_l_lo_el_n_p1, r_n, wk_a(r_n)
C? endif