-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmc64ad_dist.c
2654 lines (2543 loc) · 68.6 KB
/
mc64ad_dist.c
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
/* mc64ad.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "superlu_ddefs.h"
#define abs(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) < (b)) ? (a) : (b)
/* Table of constant values */
static int_t c__1 = 1;
static int_t c__2 = 2;
/*! @file
* \brief Permute large entries to the main diagonal
*/
/* CCCC COPYRIGHT (c) 1999 Council for the Central Laboratory of the */
/* CCCC Research Councils. All rights reserved. */
/* CCCC PACKAGE MC64A/AD */
/* CCCC AUTHORS Iain Duff ([email protected]) and Jacko Koster ([email protected]) */
/* CCCC LAST UPDATE 20/09/99 */
/* CCCC */
/* *** Conditions on external use *** */
/* The user shall acknowledge the contribution of this */
/* package in any publication of material dependent upon the use of */
/* the package. The user shall use reasonable endeavours to notify */
/* the authors of the package of this publication. */
/* The user can modify this code but, at no time */
/* shall the right or title to all or any part of this package pass */
/* to the user. The user shall make available free of charge */
/* to the authors for any purpose all information relating to any */
/* alteration or addition made to this package for the purposes of */
/* extending the capabilities or enhancing the performance of this */
/* package. */
/* The user shall not pass this code directly to a third party without the */
/* express prior consent of the authors. Users wanting to licence their */
/* own copy of these routines should send email to [email protected] */
/* None of the comments from the Copyright notice up to and including this */
/* one shall be removed or altered in any way. */
/* ********************************************************************** */
/* </pre>
*/
/* Subroutine */ int_t mc64id_dist(int_t *icntl)
{
int_t i__;
/* *** Copyright (c) 1999 Council for the Central Laboratory of the */
/* Research Councils *** */
/* *** Although every effort has been made to ensure robustness and *** */
/* *** reliability of the subroutines in this MC64 suite, we *** */
/* *** disclaim any liability arising through the use or misuse of *** */
/* *** any of the subroutines. *** */
/* *** Any problems? Contact ... */
/* Iain Duff ([email protected]) or Jacko Koster ([email protected]) *** */
/* Purpose */
/* ======= */
/* The components of the array ICNTL control the action of MC64A/AD. */
/* Default values for these are set in this subroutine. */
/* Parameters */
/* ========== */
/* Local variables */
/* ICNTL(1) has default value 6. */
/* It is the output stream for error messages. If it */
/* is negative, these messages will be suppressed. */
/* ICNTL(2) has default value 6. */
/* It is the output stream for warning messages. */
/* If it is negative, these messages are suppressed. */
/* ICNTL(3) has default value -1. */
/* It is the output stream for monitoring printing. */
/* If it is negative, these messages are suppressed. */
/* ICNTL(4) has default value 0. */
/* If left at the default value, the incoming data is checked for */
/* out-of-range indices and duplicates. Setting ICNTL(4) to any */
/* other will avoid the checks but is likely to cause problems */
/* later if out-of-range indices or duplicates are present. */
/* The user should only set ICNTL(4) non-zero, if the data is */
/* known to avoid these problems. */
/* ICNTL(5) to ICNTL(10) are not used by MC64A/AD but are set to */
/* zero in this routine. */
/* Initialization of the ICNTL array. */
/* Parameter adjustments */
--icntl;
/* Function Body */
icntl[1] = 6;
icntl[2] = 6;
icntl[3] = -1;
for (i__ = 4; i__ <= 10; ++i__) {
icntl[i__] = 0;
/* L10: */
}
return 0;
} /* mc64id_ */
/* ********************************************************************** */
/* Subroutine */ int_t mc64ad_dist(int_t *job, int_t *n, int_t *ne, int_t *
ip, int_t *irn, double *a, int_t *num, int_t *cperm,
int_t *liw, int_t *iw, int_t *ldw, double *dw, int_t *
icntl, int_t *info)
{
/* System generated locals */
int_t i__1, i__2;
double d__1, d__2;
/* Builtin functions */
double log(double);
/* Local variables */
int_t i__, j, k;
double fact, rinf;
extern /* Subroutine */ int_t mc21ad_dist(int_t *, int_t *, int_t *,
int_t *, int_t *, int_t *, int_t *, int_t *),
mc64bd_dist(int_t *, int_t *, int_t *, int_t *, double *, int_t
*, int_t *, int_t *, int_t *, int_t *, int_t *, double *),
mc64rd_dist(int_t *, int_t *, int_t *, int_t *, double *),
mc64sd_dist(int_t *, int_t *, int_t *, int_t *
, double *, int_t *, int_t *, int_t *, int_t *,
int_t *, int_t *, int_t *, int_t *, int_t *),
mc64wd_dist(int_t *, int_t *, int_t *, int_t *, double *, int_t
*, int_t *, int_t *, int_t *, int_t *, int_t *, int_t
*, double *, double *);
/* *** Copyright (c) 1999 Council for the Central Laboratory of the */
/* Research Councils *** */
/* *** Although every effort has been made to ensure robustness and *** */
/* *** reliability of the subroutines in this MC64 suite, we *** */
/* *** disclaim any liability arising through the use or misuse of *** */
/* *** any of the subroutines. *** */
/* *** Any problems? Contact ... */
/* Iain Duff ([email protected]) or Jacko Koster ([email protected]) *** */
/* Purpose */
/* ======= */
/*! \brief
* <pre>
* This subroutine attempts to find a column permutation for an NxN
* sparse matrix A = {a_ij} that makes the permuted matrix have N
* entries on its diagonal.
* If the matrix is structurally nonsingular, the subroutine optionally
* returns a column permutation that maximizes the smallest element
* on the diagonal, maximizes the sum of the diagonal entries, or
* maximizes the product of the diagonal entries of the permuted matrix.
* For the latter option, the subroutine also finds scaling factors
* that may be used to scale the matrix so that the nonzero diagonal
* entries of the permuted matrix are one in absolute value and all the
* off-diagonal entries are less than or equal to one in absolute value.
* The natural logarithms of the scaling factors u(i), i=1..N, for the
* rows and v(j), j=1..N, for the columns are returned so that the
* scaled matrix B = {b_ij} has entries b_ij = a_ij * EXP(u_i + v_j).
* </pre>
*/
/* Parameters */
/* ========== */
/* JOB is an INT_T variable which must be set by the user to */
/* control the action. It is not altered by the subroutine. */
/* Possible values for JOB are: */
/* 1 Compute a column permutation of the matrix so that the */
/* permuted matrix has as many entries on its diagonal as possible. */
/* The values on the diagonal are of arbitrary size. HSL subroutine */
/* MC21A/AD is used for this. See [1]. */
/* 2 Compute a column permutation of the matrix so that the smallest */
/* value on the diagonal of the permuted matrix is maximized. */
/* See [3]. */
/* 3 Compute a column permutation of the matrix so that the smallest */
/* value on the diagonal of the permuted matrix is maximized. */
/* The algorithm differs from the one used for JOB = 2 and may */
/* have quite a different performance. See [2]. */
/* 4 Compute a column permutation of the matrix so that the sum */
/* of the diagonal entries of the permuted matrix is maximized. */
/* See [3]. */
/* 5 Compute a column permutation of the matrix so that the product */
/* of the diagonal entries of the permuted matrix is maximized */
/* and vectors to scale the matrix so that the nonzero diagonal */
/* entries of the permuted matrix are one in absolute value and */
/* all the off-diagonal entries are less than or equal to one in */
/* absolute value. See [3]. */
/* Restriction: 1 <= JOB <= 5. */
/* N is an INT_T variable which must be set by the user to the */
/* order of the matrix A. It is not altered by the subroutine. */
/* Restriction: N >= 1. */
/* NE is an INT_T variable which must be set by the user to the */
/* number of entries in the matrix. It is not altered by the */
/* subroutine. */
/* Restriction: NE >= 1. */
/* IP is an INT_T array of length N+1. */
/* IP(J), J=1..N, must be set by the user to the position in array IRN */
/* of the first row index of an entry in column J. IP(N+1) must be set */
/* to NE+1. It is not altered by the subroutine. */
/* IRN is an INT_T array of length NE. */
/* IRN(K), K=1..NE, must be set by the user to hold the row indices of */
/* the entries of the matrix. Those belonging to column J must be */
/* stored contiguously in the positions IP(J)..IP(J+1)-1. The ordering */
/* of the row indices within each column is unimportant. Repeated */
/* entries are not allowed. The array IRN is not altered by the */
/* subroutine. */
/* A is a REAL (DOUBLE PRECISION in the D-version) array of length NE. */
/* The user must set A(K), K=1..NE, to the numerical value of the */
/* entry that corresponds to IRN(K). */
/* It is not used by the subroutine when JOB = 1. */
/* It is not altered by the subroutine. */
/* NUM is an INT_T variable that need not be set by the user. */
/* On successful exit, NUM will be the number of entries on the */
/* diagonal of the permuted matrix. */
/* If NUM < N, the matrix is structurally singular. */
/* CPERM is an INT_T array of length N that need not be set by the */
/* user. On successful exit, CPERM contains the column permutation. */
/* Column CPERM(J) of the original matrix is column J in the permuted */
/* matrix, J=1..N. */
/* LIW is an INT_T variable that must be set by the user to */
/* the dimension of array IW. It is not altered by the subroutine. */
/* Restriction: */
/* JOB = 1 : LIW >= 5N */
/* JOB = 2 : LIW >= 4N */
/* JOB = 3 : LIW >= 10N + NE */
/* JOB = 4 : LIW >= 5N */
/* JOB = 5 : LIW >= 5N */
/* IW is an INT_T array of length LIW that is used for workspace. */
/* LDW is an INT_T variable that must be set by the user to the */
/* dimension of array DW. It is not altered by the subroutine. */
/* Restriction: */
/* JOB = 1 : LDW is not used */
/* JOB = 2 : LDW >= N */
/* JOB = 3 : LDW >= NE */
/* JOB = 4 : LDW >= 2N + NE */
/* JOB = 5 : LDW >= 3N + NE */
/* DW is a REAL (DOUBLE PRECISION in the D-version) array of length LDW */
/* that is used for workspace. If JOB = 5, on return, */
/* DW(i) contains u_i, i=1..N, and DW(N+j) contains v_j, j=1..N. */
/* ICNTL is an INT_T array of length 10. Its components control the */
/* output of MC64A/AD and must be set by the user before calling */
/* MC64A/AD. They are not altered by the subroutine. */
/* ICNTL(1) must be set to specify the output stream for */
/* error messages. If ICNTL(1) < 0, messages are suppressed. */
/* The default value set by MC46I/ID is 6. */
/* ICNTL(2) must be set by the user to specify the output stream for */
/* warning messages. If ICNTL(2) < 0, messages are suppressed. */
/* The default value set by MC46I/ID is 6. */
/* ICNTL(3) must be set by the user to specify the output stream for */
/* diagnostic messages. If ICNTL(3) < 0, messages are suppressed. */
/* The default value set by MC46I/ID is -1. */
/* ICNTL(4) must be set by the user to a value other than 0 to avoid */
/* checking of the input data. */
/* The default value set by MC46I/ID is 0. */
/* INFO is an INT_T array of length 10 which need not be set by the */
/* user. INFO(1) is set non-negative to indicate success. A negative */
/* value is returned if an error occurred, a positive value if a */
/* warning occurred. INFO(2) holds further information on the error. */
/* On exit from the subroutine, INFO(1) will take one of the */
/* following values: */
/* 0 : successful entry (for structurally nonsingular matrix). */
/* +1 : successful entry (for structurally singular matrix). */
/* +2 : the returned scaling factors are large and may cause */
/* overflow when used to scale the matrix. */
/* (For JOB = 5 entry only.) */
/* -1 : JOB < 1 or JOB > 5. Value of JOB held in INFO(2). */
/* -2 : N < 1. Value of N held in INFO(2). */
/* -3 : NE < 1. Value of NE held in INFO(2). */
/* -4 : the defined length LIW violates the restriction on LIW. */
/* Value of LIW required given by INFO(2). */
/* -5 : the defined length LDW violates the restriction on LDW. */
/* Value of LDW required given by INFO(2). */
/* -6 : entries are found whose row indices are out of range. INFO(2) */
/* contains the index of a column in which such an entry is found. */
/* -7 : repeated entries are found. INFO(2) contains the index of a */
/* column in which such entries are found. */
/* INFO(3) to INFO(10) are not currently used and are set to zero by */
/* the routine. */
/* References: */
/* [1] I. S. Duff, (1981), */
/* "Algorithm 575. Permutations for a zero-free diagonal", */
/* ACM Trans. Math. Software 7(3), 387-390. */
/* [2] I. S. Duff and J. Koster, (1998), */
/* "The design and use of algorithms for permuting large */
/* entries to the diagonal of sparse matrices", */
/* SIAM J. Matrix Anal. Appl., vol. 20, no. 4, pp. 889-901. */
/* [3] I. S. Duff and J. Koster, (1999), */
/* "On algorithms for permuting large entries to the diagonal */
/* of sparse matrices", */
/* Technical Report RAL-TR-1999-030, RAL, Oxfordshire, England. */
/* Local variables and parameters */
/* External routines and functions */
/* EXTERNAL FD05AD */
/* DOUBLE PRECISION FD05AD */
/* Intrinsic functions */
/* Set RINF to largest positive real number (infinity) */
/* XSL RINF = FD05AD(5) */
/* Parameter adjustments */
--cperm;
--ip;
--a;
--irn;
--iw;
--dw;
--icntl;
--info;
/* Function Body */
rinf = dmach_dist("Overflow");
/* Check value of JOB */
if (*job < 1 || *job > 5) {
info[1] = -1;
info[2] = *job;
if (icntl[1] >= 0) {
printf(" ****** Error in MC64A/AD. INFO(1) = " IFMT
" because JOB = " IFMT "\n", info[1], *job);
}
goto L99;
}
/* Check value of N */
if (*n < 1) {
info[1] = -2;
info[2] = *n;
if (icntl[1] >= 0) {
printf(" ****** Error in MC64A/AD. INFO(1) = " IFMT
" because N = " IFMT "\n", info[1], *job);
}
goto L99;
}
/* Check value of NE */
if (*ne < 1) {
info[1] = -3;
info[2] = *ne;
if (icntl[1] >= 0) {
printf(" ****** Error in MC64A/AD. INFO(1) = " IFMT
" because NE = " IFMT "\n", info[1], *job);
}
goto L99;
}
/* Check LIW */
if (*job == 1) {
k = *n * 5;
}
if (*job == 2) {
k = *n << 2;
}
if (*job == 3) {
k = *n * 10 + *ne;
}
if (*job == 4) {
k = *n * 5;
}
if (*job == 5) {
k = *n * 5;
}
if (*liw < k) {
info[1] = -4;
info[2] = k;
if (icntl[1] >= 0) {
printf(" ****** Error in MC64A/AD. INFO(1) = " IFMT
" LIW too small, must be at least " IFMT "\n", info[1], k);
}
goto L99;
}
/* Check LDW */
/* If JOB = 1, do not check */
if (*job > 1) {
if (*job == 2) {
k = *n;
}
if (*job == 3) {
k = *ne;
}
if (*job == 4) {
k = (*n << 1) + *ne;
}
if (*job == 5) {
k = *n * 3 + *ne;
}
if (*ldw < k) {
info[1] = -5;
info[2] = k;
if (icntl[1] >= 0) {
printf(" ****** Error in MC64A/AD. INFO(1) = " IFMT
" LDW too small, must be at least " IFMT "\n", info[1], k);
}
goto L99;
}
}
if (icntl[4] == 0) {
/* Check row indices. Use IW(1:N) as workspace */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
iw[i__] = 0;
/* L3: */
}
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
i__ = irn[k];
/* Check for row indices that are out of range */
if (i__ < 1 || i__ > *n) {
info[1] = -6;
info[2] = j;
if (icntl[1] >= 0) {
printf(" ****** Error in MC64A/AD. INFO(1) = " IFMT
" Column " IFMT
" contains an entry with invalid row index " IFMT "\n",
info[1], j, i__);
}
goto L99;
}
/* Check for repeated row indices within a column */
if (iw[i__] == j) {
info[1] = -7;
info[2] = j;
if (icntl[1] >= 0) {
printf(" ****** Error in MC64A/AD. INFO(1) = " IFMT
" Column " IFMT
" contains two or more entries with row index " IFMT "\n",
info[1], j, i__);
}
goto L99;
} else {
iw[i__] = j;
}
/* L4: */
}
/* L6: */
}
}
/* Print diagnostics on input */
if (icntl[3] >= 0) {
printf(" ****** Input parameters for MC64A/AD: JOB = " IFMT ","
" N = " IFMT ", NE = " IFMT "\n", *job, *n, *ne);
printf(" IP(1:N+1) = ");
for (j=1; j<=(*n+1); ++j) {
printf(IFMT, ip[j]);
if (j%8 == 0) printf("\n");
}
printf("\n IRN(1:NE) = ");
for (j=1; j<=(*ne); ++j) {
printf(IFMT, irn[j]);
if (j%8 == 0) printf("\n");
}
printf("\n");
if (*job > 1) {
printf(" A(1:NE) = ");
for (j=1; j<=(*ne); ++j) {
printf("%f14.4", a[j]);
if (j%4 == 0) printf("\n");
}
printf("\n");
}
}
/* Set components of INFO to zero */
for (i__ = 1; i__ <= 10; ++i__) {
info[i__] = 0;
/* L8: */
}
/* Compute maximum matching with MC21A/AD */
if (*job == 1) {
/* Put length of column J in IW(J) */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
iw[j] = ip[j + 1] - ip[j];
/* L10: */
}
/* IW(N+1:5N) is workspace */
#if 0
mc21ad_(n, &irn[1], ne, &ip[1], &iw[1], &cperm[1], num, &iw[*n+1]);
#else
printf(" ****** Warning from MC64A/AD. Need to link mc21ad.\n");
#endif
goto L90;
}
/* Compute bottleneck matching */
if (*job == 2) {
/* IW(1:5N), DW(1:N) are workspaces */
mc64bd_dist(n, ne, &ip[1], &irn[1], &a[1], &cperm[1], num,
&iw[1], &iw[*n + 1], &iw[(*n << 1) + 1], &iw[*n * 3 + 1],
&dw[1]);
goto L90;
}
/* Compute bottleneck matching */
if (*job == 3) {
/* Copy IRN(K) into IW(K), ABS(A(K)) into DW(K), K=1..NE */
i__1 = *ne;
for (k = 1; k <= i__1; ++k) {
iw[k] = irn[k];
dw[k] = (d__1 = a[k], abs(d__1));
/* L20: */
}
/* Sort entries in each column by decreasing value. */
mc64rd_dist(n, ne, &ip[1], &iw[1], &dw[1]);
/* IW(NE+1:NE+10N) is workspace */
mc64sd_dist(n, ne, &ip[1], &iw[1], &dw[1], &cperm[1], num,
&iw[*ne + 1], &iw[*ne + *n + 1], &iw[*ne + (*n << 1) + 1],
&iw[*ne + *n * 3 + 1], &iw[*ne + (*n << 2) + 1],
&iw[*ne + *n * 5 + 1], &iw[*ne + *n * 6 + 1]);
goto L90;
}
if (*job == 4) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
fact = 0.;
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
if ((d__1 = a[k], abs(d__1)) > fact) {
fact = (d__2 = a[k], abs(d__2));
}
/* L30: */
}
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
dw[(*n << 1) + k] = fact - (d__1 = a[k], abs(d__1));
/* L40: */
}
/* L50: */
}
/* B = DW(2N+1:2N+NE); IW(1:5N) and DW(1:2N) are workspaces */
mc64wd_dist(n, ne, &ip[1], &irn[1], &dw[(*n << 1) + 1], &cperm[1],
num, &iw[1], &iw[*n + 1], &iw[(*n << 1) + 1],
&iw[*n * 3 + 1], &iw[(*n << 2) + 1], &dw[1], &dw[*n + 1]);
goto L90;
}
if (*job == 5) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
fact = 0.;
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
dw[*n * 3 + k] = (d__1 = a[k], abs(d__1));
if (dw[*n * 3 + k] > fact) {
fact = dw[*n * 3 + k];
}
/* L60: */
}
dw[(*n << 1) + j] = fact;
if (fact != 0.) {
fact = log(fact);
} else {
fact = rinf / *n;
}
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
if (dw[*n * 3 + k] != 0.) {
dw[*n * 3 + k] = fact - log(dw[*n * 3 + k]);
} else {
dw[*n * 3 + k] = rinf / *n;
}
/* L70: */
}
/* L75: */
}
/* B = DW(3N+1:3N+NE); IW(1:5N) and DW(1:2N) are workspaces */
mc64wd_dist(n, ne, &ip[1], &irn[1], &dw[*n * 3 + 1], &cperm[1],
num, &iw[1], &iw[*n + 1], &iw[(*n << 1) + 1],
&iw[*n * 3 + 1], &iw[(*n << 2) + 1], &dw[1], &dw[*n + 1]);
if (*num == *n) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (dw[(*n << 1) + j] != 0.) {
dw[*n + j] -= log(dw[(*n << 1) + j]);
} else {
dw[*n + j] = 0.;
}
/* L80: */
}
}
/* Check size of scaling factors */
fact = log(rinf) * .5f;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (dw[j] < fact && dw[*n + j] < fact) {
goto L86;
}
info[1] = 2;
goto L90;
L86:
;
}
/* GO TO 90 */
}
L90:
if (info[1] == 0 && *num < *n) {
/* Matrix is structurally singular, return with warning */
info[1] = 1;
if (icntl[2] >= 0) {
printf(" ****** Warning from MC64A/AD. INFO(1) = " IFMT
" The matrix is structurally singular.\n", info[1]);
}
}
if (info[1] == 2) {
/* Scaling factors are large, return with warning */
if (icntl[2] >= 0) {
printf(" ****** Warning from MC64A/AD. INFO(1) = " IFMT "\n"
" Some scaling factors may be too large.\n", info[1]);
}
}
/* Print diagnostics on output */
if (icntl[3] >= 0) {
printf(" ****** Output parameters for MC64A/AD: INFO(1:2) = " IFMT IFMT "\n",
info[1], info[2]);
printf(" NUM = " IFMT, *num);
printf(" CPERM(1:N) = ");
for (j=1; j<=*n; ++j) {
printf(IFMT, cperm[j]);
if (j%8 == 0) printf("\n");
}
if (*job == 5) {
printf("\n DW(1:N) = ");
for (j=1; j<=*n; ++j) {
printf("%11.3f", dw[j]);
if (j%5 == 0) printf("\n");
}
printf("\n DW(N+1:2N) = ");
for (j=1; j<=*n; ++j) {
printf("%11.3f", dw[*n+j]);
if (j%5 == 0) printf("\n");
}
printf("\n");
}
}
/* Return from subroutine. */
L99:
return 0;
} /* mc64ad_ */
/* ********************************************************************** */
/* Subroutine */ int_t mc64bd_dist(int_t *n, int_t *ne, int_t *ip, int_t *
irn, double *a, int_t *iperm, int_t *num, int_t *jperm,
int_t *pr, int_t *q, int_t *l, double *d__)
{
/* System generated locals */
int_t i__1, i__2, i__3;
double d__1, d__2, d__3;
/* Local variables */
int_t i__, j, k;
double a0;
int_t i0, q0;
double ai, di;
int_t ii, jj, kk;
double bv;
int_t up;
double dq0;
int_t kk1, kk2;
double csp;
int_t isp, jsp, low;
double dnew;
int_t jord, qlen, idum, jdum;
double rinf;
extern /* Subroutine */ int_t mc64dd_dist(int_t *, int_t *, int_t *,
double *, int_t *, int_t *), mc64ed_dist(int_t *, int_t *,
int_t *, double *, int_t *, int_t *), mc64fd_dist(int_t *
, int_t *, int_t *, int_t *, double *, int_t *, int_t *);
/* *** Copyright (c) 1999 Council for the Central Laboratory of the */
/* Research Councils *** */
/* *** Although every effort has been made to ensure robustness and *** */
/* *** reliability of the subroutines in this MC64 suite, we *** */
/* *** disclaim any liability arising through the use or misuse of *** */
/* *** any of the subroutines. *** */
/* *** Any problems? Contact ... */
/* Iain Duff ([email protected]) or Jacko Koster ([email protected]) *** */
/* N, NE, IP, IRN are described in MC64A/AD. */
/* A is a REAL (DOUBLE PRECISION in the D-version) array of length */
/* NE. A(K), K=1..NE, must be set to the value of the entry */
/* that corresponds to IRN(K). It is not altered. */
/* IPERM is an INT_T array of length N. On exit, it contains the */
/* matching: IPERM(I) = 0 or row I is matched to column IPERM(I). */
/* NUM is INT_T variable. On exit, it contains the cardinality of the */
/* matching stored in IPERM. */
/* IW is an INT_T work array of length 4N. */
/* DW is a REAL (DOUBLE PRECISION in D-version) work array of length N. */
/* Local variables */
/* Local parameters */
/* Intrinsic functions */
/* External subroutines and/or functions */
/* EXTERNAL FD05AD,MC64DD,MC64ED,MC64FD, DMACH */
/* DOUBLE PRECISION FD05AD, DMACH */
/* Set RINF to largest positive real number */
/* XSL RINF = FD05AD(5) */
/* Parameter adjustments */
--d__;
--l;
--q;
--pr;
--jperm;
--iperm;
--ip;
--a;
--irn;
/* Function Body */
rinf = dmach_dist("Overflow");
/* Initialization */
*num = 0;
bv = rinf;
i__1 = *n;
for (k = 1; k <= i__1; ++k) {
iperm[k] = 0;
jperm[k] = 0;
pr[k] = ip[k];
d__[k] = 0.;
/* L10: */
}
/* Scan columns of matrix; */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
a0 = -1.;
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
i__ = irn[k];
ai = (d__1 = a[k], abs(d__1));
if (ai > d__[i__]) {
d__[i__] = ai;
}
if (jperm[j] != 0) {
goto L30;
}
if (ai >= bv) {
a0 = bv;
if (iperm[i__] != 0) {
goto L30;
}
jperm[j] = i__;
iperm[i__] = j;
++(*num);
} else {
if (ai <= a0) {
goto L30;
}
a0 = ai;
i0 = i__;
}
L30:
;
}
if (a0 != -1. && a0 < bv) {
bv = a0;
if (iperm[i0] != 0) {
goto L20;
}
iperm[i0] = j;
jperm[j] = i0;
++(*num);
}
L20:
;
}
/* Update BV with smallest of all the largest maximum absolute values */
/* of the rows. */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MIN */
d__1 = bv, d__2 = d__[i__];
bv = min(d__1,d__2);
/* L25: */
}
if (*num == *n) {
goto L1000;
}
/* Rescan unassigned columns; improve initial assignment */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (jperm[j] != 0) {
goto L95;
}
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
i__ = irn[k];
ai = (d__1 = a[k], abs(d__1));
if (ai < bv) {
goto L50;
}
if (iperm[i__] == 0) {
goto L90;
}
jj = iperm[i__];
kk1 = pr[jj];
kk2 = ip[jj + 1] - 1;
if (kk1 > kk2) {
goto L50;
}
i__3 = kk2;
for (kk = kk1; kk <= i__3; ++kk) {
ii = irn[kk];
if (iperm[ii] != 0) {
goto L70;
}
if ((d__1 = a[kk], abs(d__1)) >= bv) {
goto L80;
}
L70:
;
}
pr[jj] = kk2 + 1;
L50:
;
}
goto L95;
L80:
jperm[jj] = ii;
iperm[ii] = jj;
pr[jj] = kk + 1;
L90:
++(*num);
jperm[j] = i__;
iperm[i__] = j;
pr[j] = k + 1;
L95:
;
}
if (*num == *n) {
goto L1000;
}
/* Prepare for main loop */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
d__[i__] = -1.;
l[i__] = 0;
/* L99: */
}
/* Main loop ... each pass round this loop is similar to Dijkstra's */
/* algorithm for solving the single source shortest path problem */
i__1 = *n;
for (jord = 1; jord <= i__1; ++jord) {
if (jperm[jord] != 0) {
goto L100;
}
qlen = 0;
low = *n + 1;
up = *n + 1;
/* CSP is cost of shortest path to any unassigned row */
/* ISP is matrix position of unassigned row element in shortest path */
/* JSP is column index of unassigned row element in shortest path */
csp = -1.;
/* Build shortest path tree starting from unassigned column JORD */
j = jord;
pr[j] = -1;
/* Scan column J */
i__2 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__2; ++k) {
i__ = irn[k];
dnew = (d__1 = a[k], abs(d__1));
if (csp >= dnew) {
goto L115;
}
if (iperm[i__] == 0) {
/* Row I is unassigned; update shortest path info */
csp = dnew;
isp = i__;
jsp = j;
if (csp >= bv) {
goto L160;
}
} else {
d__[i__] = dnew;
if (dnew >= bv) {
/* Add row I to Q2 */
--low;
q[low] = i__;
} else {
/* Add row I to Q, and push it */
++qlen;
l[i__] = qlen;
mc64dd_dist(&i__, n, &q[1], &d__[1], &l[1], &c__1);
}
jj = iperm[i__];
pr[jj] = j;
}
L115:
;
}
i__2 = *num;
for (jdum = 1; jdum <= i__2; ++jdum) {
/* If Q2 is empty, extract new rows from Q */
if (low == up) {
if (qlen == 0) {
goto L160;
}
i__ = q[1];
if (csp >= d__[i__]) {
goto L160;
}
bv = d__[i__];
i__3 = *n;
for (idum = 1; idum <= i__3; ++idum) {
mc64ed_dist(&qlen, n, &q[1], &d__[1], &l[1], &c__1);
l[i__] = 0;
--low;
q[low] = i__;
if (qlen == 0) {
goto L153;
}
i__ = q[1];
if (d__[i__] != bv) {
goto L153;
}
/* L152: */
}
/* End of dummy loop; this point is never reached */
}
/* Move row Q0 */
L153:
--up;
q0 = q[up];
dq0 = d__[q0];
l[q0] = up;
/* Scan column that matches with row Q0 */
j = iperm[q0];
i__3 = ip[j + 1] - 1;
for (k = ip[j]; k <= i__3; ++k) {
i__ = irn[k];
/* Update D(I) */
if (l[i__] >= up) {
goto L155;
}
/* Computing MIN */
d__2 = dq0, d__3 = (d__1 = a[k], abs(d__1));
dnew = min(d__2,d__3);
if (csp >= dnew) {
goto L155;
}
if (iperm[i__] == 0) {
/* Row I is unassigned; update shortest path info */
csp = dnew;
isp = i__;
jsp = j;
if (csp >= bv) {
goto L160;
}
} else {
di = d__[i__];
if (di >= bv || di >= dnew) {
goto L155;
}
d__[i__] = dnew;
if (dnew >= bv) {
/* Delete row I from Q (if necessary); add row I to Q2 */
if (di != -1.) {
mc64fd_dist(&l[i__], &qlen, n, &q[1], &d__[1], &l[1],
&c__1);
}
l[i__] = 0;
--low;
q[low] = i__;
} else {
/* Add row I to Q (if necessary); push row I up Q */
if (di == -1.) {
++qlen;
l[i__] = qlen;
}
mc64dd_dist(&i__, n, &q[1], &d__[1], &l[1], &c__1);
}
/* Update tree */
jj = iperm[i__];