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
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
|
SUBROUTINE TWODQ(F,N,X,Y,TOL,ICLOSE,MAXTRI,MEVALS,RESULT,
* ERROR,NU,ND,NEVALS,IFLAG,DATA,IWORK)
C***BEGIN PROLOGUE TWODQ
C***DATE WRITTEN 840518 (YYMMDD)
C***REVISION DATE 840518 (YYMMDD)
C***CATEGORY NO. D1I
C***KEYWORDS QUADRATURE,TWO DIMENSIONAL,ADAPTIVE,CUBATURE
C***AUTHOR KAHANER,D.K.,N.B.S.
C RECHARD,O.W.,N.B.S.
C BARNHILL,ROBERT,UNIV. OF UTAH
C***PURPOSE To compute the two-dimensional integral of a function
C F over a region consisting of N triangles.
C***DESCRIPTION
C
C This subroutine computes the two-dimensional integral of a
C function F over a region consisting of N triangles.
C A total error estimate is obtained and compared with a
C tolerance - TOL - that is provided as input to the subroutine.
C The error tolerance is treated as either relative or absolute
C depending on the input value of IFLAG. A 'Local Quadrature
C Module' is applied to each input triangle and estimates of the
C total integral and the total error are computed. The local
C quadrature module is either subroutine LQM0 or subroutine
C LQM1 and the choice between them is determined by the
C value of the input variable ICLOSE.
C
C If the total error estimate exceeds the tolerance, the triangle
C with the largest absolute error is divided into two triangles
C by a median to its longest side. The local quadrature module
C is then applied to each of the subtriangles to obtain new
C estimates of the integral and the error. This process is
C repeated until either (1) the error tolerance is satisfied,
C (2) the number of triangles generated exceeds the input
C parameter MAXTRI, (3) the number of integrand evaluations
C exceeds the input parameter MEVALS, or (4) the subroutine
C senses that roundoff error is beginning to contaminate
C the result.
C
C The user must specify MAXTRI, the maximum number of triangles
C in the final triangulation of the region, and provide two
C storage arrays - DATA and IWORK - whose sizes are at least
C 9*MAXTRI and 2*MAXTRI respectively. The user must also
C specify MEVALS, the maximum number of function evaluations
C to be allowed. This number will be effective in limiting
C the computation only if it is less than 94*MAXTRI when LQM1
C is specified or 56*MAXTRI when LQM0 is specified.
C
C After the subroutine has returned to the calling program
C with output values, it can be called again with a smaller
C value of TOL and/or a different value of MEVALS. The tolerance
C can also be changed from relative to absolute
C or vice-versa by changing IFLAG. Unless
C the parameters NU and ND are reset to zero the subroutine
C will restart with the final set of triangles and output
C values from the previous call.
C
C ARGUMENTS:
C
C F function subprogram defining the integrand F(u,v);
C the actual name for F needs to be declared EXTERNAL
C by the calling program.
C
C N the number of input triangles.
C
C X a 3 by N array containing the abscissae of the vertices
C of the N triangles.
C
C Y a 3 by N array containing the ordinates of the vertices
C of the N triangles
C
C TOL the desired bound on the error. If IFLAG=0 on input,
C TOL is interpreted as a bound on the relative error;
C if IFLAG=1, the bound is on the absolute error.
C
C ICLOSE an integer parameter that determines the selection
C of LQM0 or LQM1. If ICLOSE=1 then LQM1 is used.
C Any other value of ICLOSE causes LQM0 to be used.
C LQM0 uses function values only at interior points of
C the triangle. LQM1 is usually more accurate than LQM0
C but involves evaluating the integrand at more points
C including some on the boundary of the triangle. It
C will usually be better to use LQM1 unless the integrand
C has singularities on the boundary of the triangle.
C
C MAXTRI The maximum number of triangles that are allowed
C to be generated by the computation.
C
C MEVALS The maximum number of function evaluations allowed.
C
C RESULT output of the estimate of the integral.
C
C ERROR output of the estimate of the absolute value of the
C total error.
C
C NU an integer variable used for both input and output. Must
C be set to 0 on first call of the subroutine. Subsequent
C calls to restart the subroutine should use the previous
C output value.
C
C ND an integer variable used for both input and output. Must
C be set to 0 on first call of the subroutine. Subsequent
C calls to restart the subroutine should use the previous
C output value.
C
C NEVALS The actual number of function evaluations performed.
C
C IFLAG on input:
C IFLAG=0 means TOL is a bound on the relative error;
C IFLAG=1 means TOL is a bound on the absolute error;
C any other input value for IFLAG causes the subroutine
C to return immediately with IFLAG set equal to 9.
C
C on output:
C IFLAG=0 means normal termination;
C IFLAG=1 means termination for lack of space to divide
C another triangle;
C IFLAG=2 means termination because of roundoff noise
C IFLAG=3 means termination with relative error <=
C 5.0* machine epsilon;
C IFLAG=4 means termination because the number of function
C evaluations has exceeded MEVALS.
C IFLAG=9 means termination because of error in input flag
C
C DATA a one dimensional array of length >= 9*MAXTRI
C passed to the subroutine by the calling program. It is
C used by the subroutine to store information
C about triangles used in the quadrature.
C
C IWORK a one dimensional integer array of length >= 2*MAXTRI
C passed to the subroutine by the calling program.
C It is used by the subroutine to store pointers
C to the information in the DATA array.
C
C
C The information for each triangle is contained in a nine word
C record consisting of the error estimate, the estimate of the
C integral, the coordinates of the three vertices, and the area.
C These records are stored in the DATA array
C that is passed to the subroutine. The storage is organized
C into two heaps of length NU and ND respectively. The first heap
C contains those triangles for which the error exceeds
C epsabs*a/ATOT where epsabs is a bound on the absolute error
C derived from the input tolerance (which may refer to relative
C or absolute error), a is the area of the triangle, and ATOT
C is the total area of all triangles. The second heap contains
C those triangles for which the error is less than or equal to
C epsabs*a/ATOT. At the top of each heap is the triangle with
C the largest absolute error.
C
C Pointers into the heaps are contained in the array IWORK.
C Pointers to the first heap are contained
C between IWORK(1) and IWORK(NU). Pointers to the second
C heap are contained between IWORK(MAXTRI+1) and
C IWORK(MAXTRI+ND). The user thus has access to the records
C stored in the DATA array through the pointers in IWORK.
C For example, the following two DO loops will print out
C the records for each triangle in the two heaps:
C
C DO 10 I=1,NU
C PRINT*,(DATA(IWORK(I)+J),J=0,8)
C 10 CONTINUE
C DO 20 I=1,ND
C PRINT*,(DATA(IWORK(MAXTRI+I)+J),J=0,8
C 20 CONTINUE
C
C When the total number of triangles is equal to
C MAXTRI, the program attempts to remove a triangle from the
C bottom of the second heap and continue. If the second heap
C is empty, the program returns with the current estimates of
C the integral and the error and with IFLAG set equal to 1.
C Note that in this case the actual number of triangles
C processed may exceed MAXTRI and the triangles stored in
C the DATA array may not constitute a complete triangulation
C of the region.
C
C The following sample program will calculate the integral of
C cos(x+y) over the square (0.,0.),(1.,0.),(1.,1.),(0.,1.) and
C print out the values of the estimated integral, the estimated
C error, the number of function evaluations, and IFLAG.
C
C Double precision X(3,2),Y(3,2),DATA(450),RES,ERR
C INTEGER IWORK(100),NU,ND,NEVALS,IFLAG
C EXTERNAL F
C X(1,1)=0.
C Y(1,1)=0.
C X(2,1)=1.
C Y(2,1)=0.
C X(3,1)=1.
C Y(3,1)=1.
C X(1,2)=0.
C Y(1,2)=0.
C X(2,2)=1.
C Y(2,2)=1.
C X(3,2)=0.
C Y(3,2)=1.
C NU=0
C ND=0
C IFLAG=1
C CALL TWODQ(F,2,X,Y,1.D-04,1,50,4000,RES,ERR,NU,ND,
C * NEVALS,IFLAG,DATA,IWORK)
C PRINT*,RES,ERR,NEVALS,IFLAG
C END
C DOUBLE PRECISION FUNCTION F(X,Y)
C DOUBLE PRECISION X,Y
C F=COS(X+Y)
C RETURN
C END
C
C***REFERENCES (NONE)
C
C***ROUTINES CALLED HINITD,HINITU,HPACC,HPDEL,HPINS,LQM0,LQM1,
C TRIDV,DLAMCH
C***END PROLOGUE TWODQ
integer n,iflag,nevals,iclose,nu,nd,mevals,iwork(*),maxtri
double precision f,x(3,n),y(3,n),data(*),tol,result,error
integer rndcnt
logical full
double precision a,r,e,u(3),v(3),node(9),node1(9),node2(9),
* epsabs,EMACH,DLAMCH,ATOT,fadd,newres,newerr
external f,GREATR
common/iertwo/iero
save ATOT
EMACH=DLAMCH('p')
c
c If heaps are empty, apply LQM to each input triangle and
c place all of the data on the second heap.
c
if((nu+nd).eq.0) then
call HINITU(maxtri,9,nu,iwork)
call HINITD(maxtri,9,nd,iwork(maxtri+1))
ATOT=0.0
result=0.0
error=0.0
rndcnt=0
nevals=0
do 20 i=1,n
do 10 j=1,3
u(j)=x(j,i)
v(j)=y(j,i)
10 continue
a=0.5*abs(u(1)*v(2)+u(2)*v(3)+u(3)*v(1)
1 -u(1)*v(3)-u(2)*v(1)-u(3)*v(2))
ATOT=ATOT+a
if(iclose.eq.1) then
call LQM1(f,u,v,r,e)
if(iero.ne.0) return
nevals=nevals+47
else
call LQM0(f,u,v,r,e)
if(iero.ne.0) return
nevals=nevals+28
end if
result=result+r
error=error+e
node(1)=e
node(2)=r
node(3)=x(1,i)
node(4)=y(1,i)
node(5)=x(2,i)
node(6)=y(2,i)
node(7)=x(3,i)
node(8)=y(3,i)
node(9)=a
call HPINS(maxtri,9,data,nd,iwork(maxtri+1),node,GREATR)
20 continue
end if
c
c Check that input tolerance is consistent with
c machine epsilon.
c
if(iflag.eq.0) then
if(tol.le.5.0*EMACH) then
tol=5.0*EMACH
fadd=3
else
fadd=0
end if
epsabs=tol*abs(result)
else if(iflag.eq.1) then
if(tol.le.5.0*EMACH*abs(result)) then
epsabs=5.0*EMACH*abs(result)
else
fadd=0
epsabs=tol
end if
else
iflag=9
return
end if
c
c Adjust the second heap on the basis of the current
c value of epsabs.
c
2 if(nd.eq.0) go to 40
j=nd
3 if(j.eq.0) go to 40
call HPACC(maxtri,9,data,nd,iwork(maxtri+1),node,j)
if(node(1).gt.epsabs*node(9)/ATOT) then
call HPINS(maxtri,9,data,nu,iwork,node,GREATR)
call HPDEL(maxtri,9,data,nd,iwork(maxtri+1),GREATR,j)
if(j.gt.nd) j=j-1
else
j=j-1
endif
go to 3
c
c Beginning of main loop from here to end
c
40 if(nevals.ge.mevals) then
iflag=4
return
end if
if(error.le.epsabs) then
if(iflag.eq.0) then
if(error.le.abs(result)*tol) then
iflag=fadd
return
else
epsabs=abs(result)*tol
go to 2
end if
else
if(error.le.tol) then
iflag=0
return
else if(error.le.5.0*EMACH*abs(result)) then
iflag=3
return
else
epsabs=5.0*EMACH*abs(result)
go to 2
end if
end if
end if
c
c If there are too many triangles and second heap
c is not empty remove bottom triangle from second
c heap. If second heap is empty return with iflag
c set to 1 or 4.
c
if((nu+nd).ge.maxtri) then
full=.true.
if(nd.gt.0) then
iwork(nu+1)=iwork(maxtri+nd)
nd=nd-1
else
iflag=1
return
end if
else
full=.false.
end if
c
c Find triangle with largest error, divide it in
c two, and apply LQM to each half.
c
if(nd.eq.0) then
call HPACC(maxtri,9,data,nu,iwork,node,1)
call HPDEL(maxtri,9,data,nu,iwork,GREATR,1)
else if(nu.eq.0) then
call HPACC(maxtri,9,data,nd,iwork(maxtri+1),node,1)
call HPDEL(maxtri,9,data,nd,iwork(maxtri+1),GREATR,1)
else if(data(iwork(1)).ge.data(iwork(maxtri+1))) then
if(full) iwork(maxtri+nd+2)=iwork(nu)
call HPACC(maxtri,9,data,nu,iwork,node,1)
call HPDEL(maxtri,9,data,nu,iwork,GREATR,1)
else
if(full) iwork(nu+2)=iwork(maxtri+nd)
call HPACC(maxtri,9,data,nd,iwork(maxtri+1),node,1)
call HPDEL(maxtri,9,data,nd,iwork(maxtri+1),GREATR,1)
end if
call tridv(node,node1,node2,0.5d0,1)
do 60 j=1,3
u(j)=node1(2*j+1)
v(j)=node1(2*j+2)
60 continue
if(iclose.eq.1) then
call LQM1(f,u,v,node1(2),node1(1))
if(iero.ne.0) return
nevals=nevals+47
else
call LQM0(f,u,v,node1(2),node1(1))
if(iero.ne.0) return
nevals=nevals+28
end if
do 70 j=1,3
u(j)=node2(2*j+1)
v(j)=node2(2*j+2)
70 continue
if(iclose.eq.1) then
call LQM1(f,u,v,node2(2),node2(1))
if(iero.ne.0) return
nevals=nevals+47
else
call LQM0(f,u,v,node2(2),node2(1))
if(iero.ne.0) return
nevals=nevals+28
end if
newerr=node1(1)+node2(1)
newres=node1(2)+node2(2)
if(newerr.gt.0.99*node(1)) then
if(abs(node(2)-newres).le.1.D-04*abs(newres)) rndcnt=rndcnt+1
end if
result=result-node(2)+newres
error=error-node(1)+newerr
if(node1(1).gt.node1(9)*epsabs/ATOT) then
call HPINS(maxtri,9,data,nu,iwork,node1,GREATR)
else
call HPINS(maxtri,9,data,nd,iwork(maxtri+1),node1,GREATR)
end if
if(node2(1).gt.node2(9)*epsabs/ATOT) then
call HPINS(maxtri,9,data,nu,iwork,node2,GREATR)
else
call HPINS(maxtri,9,data,nd,iwork(maxtri+1),node2,GREATR)
end if
if(rndcnt.ge.20) then
iflag=2
return
end if
if(iflag.eq.0) then
if(epsabs.lt.0.5*tol*abs(result)) then
epsabs=tol*abs(result)
j=nu
5 if(j.eq.0) go to 40
call HPACC(maxtri,9,data,nu,iwork,node,j)
if(node(1).le.epsabs*node(9)/ATOT) then
call HPINS(maxtri,9,data,nd,iwork(maxtri+1),node,GREATR)
call HPDEL(maxtri,9,data,nu,iwork,GREATR,j)
if(j.gt.nu) j=j-1
else
j=j-1
end if
go to 5
end if
end if
go to 40
end
LOGICAL FUNCTION GREATR(A,B,NWDS)
INTEGER NWDS
DOUBLE PRECISION A(NWDS), B(NWDS)
GREATR= A(1) .GT. B(1)
RETURN
END
SUBROUTINE HINITD(NMAX,NWDS,N,T)
C PURPOSE
C THIS ROUTINE INITIALIZES THE HEAP PROGRAMS WITH T(NMAX)
C POINTING TO THE TOP OF THE HEAP.
C IT IS CALLED ONCE AT THE START OF EACH NEW CALCULATION.
C INPUT
C NMAX=MAXIMUM NUMBER OF NODES ALLOWED BY USER
C NWDS=NUMBER OF WORDS PER NODE
C OUTPUT
C N=CURRENT NUMBER OF NODES IN HEAP = 0.
C T=INTEGER ARRAY OF POINTERS TO POTENTIAL HEAP NODES.
C
INTEGER T(1)
DO 1 I=1,NMAX
1 T(I)=(NMAX-I)*NWDS+1
N=0
RETURN
END
SUBROUTINE HINITU(NMAX,NWDS,N,T)
C H E A P PACKAGE
C A COLLECTION OF PROGRAMS WHICH MAINTAIN A HEAP DATA
C STRUCTURE. BY CALLING THESE SUBROUTINES IT IS POSSIBLE TO
C INSERT, DELETE, AND ACCESS AN EXISTING HEAP OR TO BUILD A
C NEW HEAP FROM AN UNORDERED COLLECTION OF NODES. THE HEAP
C FUNCTION IS AN ARGUMENT TO THE SUBROUTINES ALLOWING VERY
C GENERAL ORGANIZATIONS.
C THE USER MUST DECIDE ON THE MAXIMUM NUMBER OF NODES
C ALLOWED AND DIMENSION THE REAL ARRAY DATA AND THE INTEGER
C ARRAY T USED INTERNALLY BY THE PACKAGE. THESE VARIABLES ARE
C THEN PASSED THROUGH THE CALL SEQUENCE BETWEEN THE HEAP
C PROGRAMS BUT ARE NOT IN GENERAL ACCESSED BY THE USER. HE
C MUST ALSO PROVIDE A HEAP FUNCTION WHOSE NAME MUST BE INCLUD-
C ED IN AN EXTERNAL STATEMENT IN THE USER PROGRAM WHICH CALLS
C THE HEAP SUBROUTINES. TWO SIMPLE HEAP FUNCTIONS ARE
C PROVIDED WITH THE PACKAGE.
C
C
C PURPOSE
C THIS ROUTINE INITIALIZES THE HEAP PROGRAMS WITH T(1)
C POINTING TO THE TOP OF THE HEAP.
C IT IS CALLED ONCE AT THE START OF EACH NEW CALCULATION
C INPUT
C NMAX = MAXIMUM NUMBER OF NODES ALLOWED BY USER.
C NWDS = NUMBER OF WORDS PER NODE
C OUTPUT
C N = CURRENT NUMBER OF NODES IN HEAP = 0.
C T = INTEGER ARRAY OF POINTERS TO POTENTIAL HEAP NODES.
C
INTEGER T(1)
DO 1 I = 1, NMAX
1 T(I)=(I-1)*NWDS+1
N = 0
RETURN
END
SUBROUTINE HPACC(NMAX,NWDS,DATA,N,T,XNODE,K)
C
C PURPOSE
C TO ACCESS THE K-TH NODE OF THE HEAP,
C 1 .LE. K .LE. N .LE. NMAX
C INPUT
C NMAX = MAXIMUM NUMBER OF NODES ALLOWED BY USER.
C DATA = WORK AREA FOR STORING NODES.
C N = CURRENT NUMBER OF NODES IN THE HEAP.
C T = INTEGER ARRAY OF POINTERS TO HEAP NODES.
C XNODE = A REAL ARRAY, NWDS WORDS LONG, IN WHICH NODAL IN-
C FORMATION WILL BE INSERTED.
C K = THE INDEX OF THE NODE TO BE FOUND AND INSERTED INTO
C XNODE.
C
C OUTPUT
C XNODE = A REAL ARRAY. CONTAINS IN XNODE(1),...,XNODE(NWDS)
C THE ELEMENTS OF THE K-TH NODE.
C
DOUBLE PRECISION DATA(1), XNODE(1)
INTEGER T(1)
IF (K .LT. 1 .OR. K .GT. N .OR. N .GT. NMAX) RETURN
J=T(K)-1
DO 1 I=1,NWDS
IPJ=I+J
1 XNODE(I)=DATA(IPJ)
RETURN
END
SUBROUTINE HPDEL(NMAX,NWDS,DATA,N,T,HPFUN,K)
C
C PURPOSE
C DELETE K-TH ELEMENT OF HEAP. RESULTING TREE IS REHEAPED.
C INPUT
C NMAX = MAXIMUN NUMBER OF NODES ALLOWED BY USER.
C NWDS = NUMBER OF WORDS PER NODE.
C DATA = WORK AREA IN WHICH THE NODES ARE STORED.
C N = CURRENT NUMBER OF NODES.
C T = INTEGER ARRAY OF POINTERS TO NODES.
C HPFUN = NAME OF USER WRITTEN FUNCTION TO DETERMINE TOP NODE.
C K = INDEX OF NODE TO BE DELETED
C OUTPUT
C N = UPDATED NUMBER OF NODES.
C T = UPDATED INTEGER POINTER ARRAY TO NODES.
C
EXTERNAL HPFUN
LOGICAL HPFUN
DOUBLE PRECISION DATA(1)
INTEGER T(1)
IF(N .EQ. 0) RETURN
IF(K.EQ.N) THEN
N=N-1
RETURN
END IF
KDEL=K
JUNK=T(KDEL)
T(KDEL)=T(N)
T(N)=JUNK
N=N-1
10 IF(KDEL.EQ.1) THEN
CALL HPGRO(NMAX,NWDS,DATA,N,T,HPFUN,KDEL)
RETURN
ELSE
KHALVE=KDEL/2
IL=T(KHALVE)
IR=T(KDEL)
IF(HPFUN(DATA(IL),DATA(IR),NWDS)) THEN
CALL HPGRO(NMAX,NWDS,DATA,N,T,HPFUN,KDEL)
RETURN
ELSE
T(KHALVE)=IR
T(KDEL)=IL
KDEL=KHALVE
END IF
END IF
GO TO 10
END
SUBROUTINE HPGRO(NMAX,NWDS,DATA,N,T,HPFUN,I)
C
C PURPOSE
C FORMS A HEAP OUT OF A TREE. USED PRIVATELY BY HPBLD.
C THE TOP OF THE TREE IS STORED IN LOCATION T(I).
C FIRST SON IS IN LOCATION T(2I), NEXT SON
C IS IN LOCATION T(2I+1).
C THIS PROGRAM ASSUMES EACH BRANCH OF THE TREE IS A HEAP.
C
INTEGER T(1)
DOUBLE PRECISION DATA(1)
LOGICAL HPFUN
IF(N .GT. NMAX) RETURN
C
K=I
1 J=2*K
C
C TEST IF ELEMENT IN J TH POSITION IS A LEAF.
C
IF( J .GT. N ) RETURN
C
C IF THERE IS MORE THAN ONE SON, FIND WHICH SON IS SMALLEST.
C
IF( J .EQ. N ) GO TO 2
IR=T(J)
IL=T(J+1)
IF(HPFUN(DATA(IL),DATA(IR),NWDS)) J=J+1
C
C IF A SON IS LARGER THAN FATHER, INTERCHANGE
C THIS DESTROYS HEAP PROPERTY, SO MUST RE-HEAP REMAINING
C ELEMENTS
C
2 CONTINUE
IL=T(K)
IR=T(J)
IF(HPFUN(DATA(IL),DATA(IR),NWDS)) RETURN
ITEMP=T(J)
T(J)=T(K)
T(K)=ITEMP
K=J
GO TO 1
END
SUBROUTINE HPINS(NMAX,NWDS,DATA,N,T,XNODE,HPFUN)
C
C PURPOSE
C THIS ROUTINE INSERTS A NODE INTO AN ALREADY EXISTING HEAP.
C THE RESULTING TREE IS RE-HEAPED.
C
C INPUT
C NMAX = MAXIMUM NUMBER OF NODES ALLOWED BY USER.
C NWDS = NUMBER OF WORDS PER NODE.
C DATA = WORK AREA FOR STORING NODES.
C N = CURRENT NUMBER OF NODES IN THE TREE.
C T = INTEGER ARRAY OF POINTERS TO HEAP NODES.
C XNODE = A REAL ARRAY, NWDS WORDS LONG, WHICH
C CONTAINS THE NODAL INFORMATION TO BE INSERTED.
C HPFUN = NAME OF USER WRITTEN FUNCTION TO DETERMINE
C THE TOP NODE.
C OUTPUT
C DATA = WORK AREA WITH NEW NODE INSERTED.
C N = UPDATED NUMBER OF NODES.
C T = UPDATED INTEGER POINTER ARRAY.
C
DOUBLE PRECISION XNODE(1),DATA(1)
INTEGER T(1)
LOGICAL HPFUN
EXTERNAL HPFUN
IF(N .EQ. NMAX) RETURN
N=N+1
J= T(N)-1
DO 1 I= 1,NWDS
IPJ=I+J
1 DATA(IPJ) = XNODE(I)
J=N
2 CONTINUE
IF(J .EQ. 1) RETURN
JR=T(J)
J2=J/2
JL=T(J2)
IF(HPFUN(DATA(JL),DATA(JR),NWDS)) RETURN
T(J2)=T(J)
T(J)=JL
J=J2
GO TO 2
END
SUBROUTINE LQM0(F,U,V,RES8,EST)
C
C
C
C PURPOSE
C TO COMPUTE - IF = INTEGRAL OF F OVER THE TRIANGLE
C WITH VERTICES (U(1),V(1)),(U(2),V(2)),(U(3),V(3)), AND
C ESTIMATE THE ERROR,
C - INTEGRAL OF ABS(F) OVER THIS TRIANGLE
C
C CALLING SEQUENCE
C CALL LQM0(F,U,V,RES11,EST)
C PARAMETERS
C F - FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C F(X,Y); THE ACTUAL NAME FOR F NEEDS TO BE
C DECLARED E X T E R N A L IN THE CALLING
C PROGRAM
C U(1),U(2),U(3)- ABSCISSAE OF VERTICES
C V(1),V(2),V(3)- ORDINATES OF VERTICES
C RES6 - APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE
C LYNESS AND JESPERSEN RULE OF DEGREE 6, USING
C 12 POINTS
C RES8 - APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE
C LYNESS AND JESPERSEN RULE OF DEGREE 8,
C USING 16 POINTS
C EST - ESTIMATE OF THE ABSOLUTE ERROR
C DRESC - APPROXIMATION TO THE INTEGRAL OF ABS(F- IF/DJ),
C OBTAINED BY THE RULE OF DEGREE 6, AND USED FOR
C THE COMPUTATION OF EST
C
C REMARKS
C DATE OF LAST UPDATE : 10 APRIL 1984 O.W. RECHARD NBS
C
C SUBROUTINES OR FUNCTIONS CALLED :
C - F (USER-SUPPLIED INTEGRAND FUNCTION)
C - DLAMCH FOR MACHINE DEPENDENT INFORMATION
C
C .....................................................................
C
DOUBLE PRECISION DJ,DF0,DRESC,EMACH,EST,F,FV,F0,U,V,
* RES6,RES8,U1,U2,U3,UFLOW,V1,V2,V3,W,W60,W80,X,Y
* ,ZETA1,ZETA2,Z1,Z2,Z3,RESAB6
EXTERNAL F
DOUBLE PRECISION DLAMCH
INTEGER J,KOUNT,L
common/iertwo/iero
C
DIMENSION FV(19),W(9),X(3),Y(3),ZETA1(9),ZETA2(9),U(3),V(3)
C
C
C FIRST HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-6
C AND DEGREE-8 FORMULA, TAKEN WITH MULTIPLICITY 3
DATA ZETA1(1),ZETA1(2),ZETA1(3),ZETA1(4),ZETA1(5),ZETA1(6),ZETA1(7
* ),ZETA1(8),ZETA1(9)/0.5014265096581342D+00,
* 0.8738219710169965D+00,0.6365024991213939D+00,
* 0.5314504984483216D-01,0.8141482341455413D-01,
* 0.8989055433659379D+00,0.6588613844964797D+00,
* 0.8394777409957211D-02,0.7284923929554041D+00/
C SECOND HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-6
C AND DEGREE-8 FORMULA, TAKEN WITH MUNLTIPLICITY 3
DATA ZETA2(1),ZETA2(2),ZETA2(3),ZETA2(4),ZETA2(5),ZETA2(6),ZETA2(7
* ),ZETA2(8),ZETA2(9)/0.2492867451709329D+00,
* 0.6308901449150177D-01,0.5314504984483216D-01,
* 0.6365024991213939D+00,0.4592925882927229D+00,
* 0.5054722831703103D-01,0.1705693077517601D+00,
* 0.7284923929554041D+00,0.8394777409957211D-02/
C WEIGHTS OF MID-POINT OF TRIANGLE IN DEGREE-6
C RESP. DEGREE-8 FORMULAE
DATA W60/0.0D+00/
DATA W80/0.1443156076777862D+00/
C WEIGHTS IN DEGREE-6 AND DEGREE-8 RULE
DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9)/
* 0.1167862757263407D+00,0.5084490637020547D-01,
* 0.8285107561839291D-01,0.8285107561839291D-01,
* 0.9509163426728497D-01,0.3245849762319813D-01,
* 0.1032173705347184D+00,0.2723031417443487D-01,
* 0.2723031417443487D-01/
C
C LIST OF MAJOR VARIABLES
C ----------------------
C DJ - AREA OF THE TRIANGLE
C DRESC - APPROXIMATION TO INTEGRAL OF
C ABS(F- IF/DJ) OVER THE TRIANGLE
C RESAB6 - APPROXIMATION TO INTEGRAL OF
C ABS(F) OVER THE TRIANGLE
C X - CARTESIAN ABSCISSAE OF THE INTEGRATION
C POINTS
C Y - CARTESIAN ORDINATES OF THE INTEGRATION
C POINTS
C FV - FUNCTION VALUES
C
C COMPUTE DEGREE-6 AND DEGREE-8 RESULTS FOR IF/DJ AND
C DEGREE-6 APPROXIMATION FOR ABS(F)
C
EMACH = DLAMCH('p')
UFLOW = DLAMCH('u')
U1=U(1)
U2=U(2)
U3=U(3)
V1=V(1)
V2=V(2)
V3=V(3)
DJ = ABS(U1*V2-U2*V1-U1*V3+V1*U3+U2*V3-V2*U3)*0.5D+00
F0 = F((U1+U2+U3)/3.0D+00,(V1+V2+V3)/3.0D+00)
if(iero.ne.0) return
RES6 = F0*W60
RESAB6 = ABS(F0)*W60
FV(1) = F0
KOUNT = 1
RES8 = F0*W80
DO 50 J=1,9
Z1 = ZETA1(J)
Z2 = ZETA2(J)
Z3 = 1.0D+00-Z1-Z2
X(1) = Z1*U1+Z2*U2+Z3*U3
Y(1) = Z1*V1+Z2*V2+Z3*V3
X(2) = Z2*U1+Z3*U2+Z1*U3
Y(2) = Z2*V1+Z3*V2+Z1*V3
X(3) = Z3*U1+Z1*U2+Z2*U3
Y(3) = Z3*V1+Z1*V2+Z2*V3
IF(J.LE.4) THEN
F0 = 0.0D+00
DF0 = 0.0D+00
DO 10 L=1,3
KOUNT = KOUNT+1
FV(KOUNT) = F(X(L),Y(L))
if(iero.ne.0) return
F0 = F0+FV(KOUNT)
DF0 = DF0+ABS(FV(KOUNT))
10 CONTINUE
RES6 = RES6+F0*W(J)
RESAB6 = RESAB6+DF0*W(J)
ELSE
F0 = F(X(1),Y(1))+F(X(2),Y(2))+F(X(3),Y(3))
if(iero.ne.0) return
RES8 = RES8+F0*W(J)
ENDIF
50 CONTINUE
C
C COMPUTE DEGREE-6 APPROXIMATION FOR THE INTEGRAL OF
C ABS(F-IF/DJ)
C
DRESC = ABS(FV(1)-RES6)*W60
KOUNT = 2
DO 60 J=1,4
DRESC = DRESC+(ABS(FV(KOUNT)-RES6)+ABS(FV(KOUNT+1)-RES6)+ABS(
* FV(KOUNT+2)-RES6))*W(J)
KOUNT = KOUNT+3
60 CONTINUE
C
C COMPUTE DEGREE-6 AND DEGREE-8 APPROXIMATIONS FOR IF,
C AND ERROR ESTIMATE
C
RES6 = RES6*DJ
RES8 = RES8*DJ
RESAB6 = RESAB6*DJ
DRESC = DRESC*DJ
EST = ABS(RES6-RES8)
IF(DRESC.NE.0.0D+00) EST = MAX(EST,DRESC*MIN(1.0D+00,(20.0D+00
* *EST/DRESC)**1.5D+00))
IF(RESAB6.GT.UFLOW) EST = MAX(EMACH*RESAB6,EST)
RETURN
END
SUBROUTINE LQM1(F,U,V,RES11,EST)
C
C
C
C PURPOSE
C TO COMPUTE - IF = INTEGRAL OF F OVER THE TRIANGLE
C WITH VERTICES (U(1),V(1)),(U(2),V(2)),(U(3),V(3)), AND
C ESTIMATE THE ERROR,
C - INTEGRAL OF ABS(F) OVER THIS TRIANGLE
C
C CALLING SEQUENCE
C CALL LQM1(F,U,V,RES11,EST)
C PARAMETERS
C F - FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C F(X,Y); THE ACTUAL NAME FOR F NEEDS TO BE
C DECLARED E X T E R N A L IN THE CALLING
C PROGRAM
C U(1),U(2),U(3)- ABSCISSAE OF VERTICES
C V(1),V(2),V(3)- ORDINATES OF VERTICES
C RES9 - APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE
C LYNESS AND JESPERSEN RULE OF DEGREE 9, USING
C 19 POINTS
C RES11 - APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE
C LYNESS AND JESPERSEN RULE OF DEGREE 11,
C USING 28 POINTS
C EST - ESTIMATE OF THE ABSOLUTE ERROR
C DRESC - APPROXIMATION TO THE INTEGRAL OF ABS(F- IF/DJ),
C OBTAINED BY THE RULE OF DEGREE 9, AND USED FOR
C THE COMPUTATION OF EST
C
C REMARKS
C DATE OF LAST UPDATE : 18 JAN 1984 D. KAHANER NBS
C
C SUBROUTINES OR FUNCTIONS CALLED :
C - F (USER-SUPPLIED INTEGRAND FUNCTION)
C - DLAMCH FOR MACHINE DEPENDENT INFORMATION
C
C .....................................................................
C
DOUBLE PRECISION DJ,DF0,DRESC,EMACH,EST,F,FV,F0,U,V,
* RES9,RES11,U1,U2,U3,UFLOW,V1,V2,V3,W,W90,W110,X,Y
* ,ZETA1,ZETA2,Z1,Z2,Z3
EXTERNAL F
DOUBLE PRECISION DLAMCH
INTEGER J,KOUNT,L
common/iertwo/iero
C
DIMENSION FV(19),W(15),X(3),Y(3),ZETA1(15),ZETA2(15),U(3),V(3)
C
C
C FIRST HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-9
C AND DEGREE-11 FORMULA, TAKEN WITH MULTIPLICITY 3
DATA ZETA1(1),ZETA1(2),ZETA1(3),ZETA1(4),ZETA1(5),ZETA1(6),ZETA1(7
* ),ZETA1(8),ZETA1(9),ZETA1(10),ZETA1(11),ZETA1(12),ZETA1(13),
* ZETA1(14),ZETA1(15)/0.2063496160252593D-01,0.1258208170141290D+
* 00,0.6235929287619356D+00,0.9105409732110941D+00,
* 0.3683841205473626D-01,0.7411985987844980D+00,
* 0.9480217181434233D+00,0.8114249947041546D+00,
* 0.1072644996557060D-01,0.5853132347709715D+00,
* 0.1221843885990187D+00,0.4484167758913055D-01,
* 0.6779376548825902D+00,0.0D+00,0.8588702812826364D+00/
C SECOND HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-9
C AND DEGREE-11 FORMULA, TAKEN WITH MUNLTIPLICITY 3
DATA ZETA2(1),ZETA2(2),ZETA2(3),ZETA2(4),ZETA2(5),ZETA2(6),ZETA2(7
* ),ZETA2(8),ZETA2(9),ZETA2(10),ZETA2(11),ZETA2(12),ZETA2(13),
* ZETA2(14),ZETA2(15)/0.4896825191987370D+00,0.4370895914929355D+
* 00,0.1882035356190322D+00,0.4472951339445297D-01,
* 0.7411985987844980D+00,0.3683841205473626D-01,
* 0.2598914092828833D-01,0.9428750264792270D-01,
* 0.4946367750172147D+00,0.2073433826145142D+00,
* 0.4389078057004907D+00,0.6779376548825902D+00,
* 0.4484167758913055D-01,0.8588702812826364D+00,0.0D+00/
C WEIGHTS OF MID-POINT OF TRIANGLE IN DEGREE-9
C RESP. DEGREE-11 FORMULAE
DATA W90/0.9713579628279610D-01/
DATA W110/0.8797730116222190D-01/
C WEIGHTS IN DEGREE-9 AND DEGREE-11 RULE
DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10),W(11),W(12
* ),W(13),W(14),W(15)/0.3133470022713983D-01,0.7782754100477543D-
* 01,0.7964773892720910D-01,0.2557767565869810D-01,
* 0.4328353937728940D-01,0.4328353937728940D-01,
* 0.8744311553736190D-02,0.3808157199393533D-01,
* 0.1885544805613125D-01,0.7215969754474100D-01,
* 0.6932913870553720D-01,0.4105631542928860D-01,
* 0.4105631542928860D-01,0.7362383783300573D-02,
* 0.7362383783300573D-02/
C
C LIST OF MAJOR VARIABLES
C ----------------------
C DJ - AREA OF THE TRIANGLE
C DRESC - APPROXIMATION TO INTEGRAL OF
C ABS(F- IF/DJ) OVER THE TRIANGLE
C RESAB9 - APPROXIMATION TO INTEGRAL OF
C ABS(F) OVER THE TRIANGLE
C X - CARTESIAN ABSCISSAE OF THE INTEGRATION
C POINTS
C Y - CARTESIAN ORDINATES OF THE INTEGRATION
C POINTS
C FV - FUNCTION VALUES
C
C COMPUTE DEGREE-9 AND DEGREE-11 RESULTS FOR IF/DJ AND
C DEGREE-9 APPROXIMATION FOR ABS(F)
C
EMACH = DLAMCH('p')
UFLOW = DLAMCH('u')
U1=U(1)
U2=U(2)
U3=U(3)
V1=V(1)
V2=V(2)
V3=V(3)
DJ = ABS(U1*V2-U2*V1-U1*V3+V1*U3+U2*V3-V2*U3)*0.5D+00
F0 = F((U1+U2+U3)/3.0D+00,(V1+V2+V3)/3.0D+00)
if(iero.ne.0) return
RES9 = F0*W90
RESAB9 = ABS(F0)*W90
FV(1) = F0
KOUNT = 1
RES11 = F0*W110
DO 50 J=1,15
Z1 = ZETA1(J)
Z2 = ZETA2(J)
Z3 = 1.0D+00-Z1-Z2
X(1) = Z1*U1+Z2*U2+Z3*U3
Y(1) = Z1*V1+Z2*V2+Z3*V3
X(2) = Z2*U1+Z3*U2+Z1*U3
Y(2) = Z2*V1+Z3*V2+Z1*V3
X(3) = Z3*U1+Z1*U2+Z2*U3
Y(3) = Z3*V1+Z1*V2+Z2*V3
IF(J.LE.6) THEN
F0 = 0.0D+00
DF0 = 0.0D+00
DO 10 L=1,3
KOUNT = KOUNT+1
FV(KOUNT) = F(X(L),Y(L))
if(iero.ne.0) return
F0 = F0+FV(KOUNT)
DF0 = DF0+ABS(FV(KOUNT))
10 CONTINUE
RES9 = RES9+F0*W(J)
RESAB9 = RESAB9+DF0*W(J)
ELSE
F0 = F(X(1),Y(1))+F(X(2),Y(2))+F(X(3),Y(3))
if(iero.ne.0) return
RES11 = RES11+F0*W(J)
ENDIF
50 CONTINUE
C
C COMPUTE DEGREE-9 APPROXIMATION FOR THE INTEGRAL OF
C ABS(F-IF/DJ)
C
DRESC = ABS(FV(1)-RES9)*W90
KOUNT = 2
DO 60 J=1,6
DRESC = DRESC+(ABS(FV(KOUNT)-RES9)+ABS(FV(KOUNT+1)-RES9)+ABS(
* FV(KOUNT+2)-RES9))*W(J)
KOUNT = KOUNT+3
60 CONTINUE
C
C COMPUTE DEGREE-9 AND DEGREE-11 APPROXIMATIONS FOR IF,
C AND ERROR ESTIMATE
C
RES9 = RES9*DJ
RES11 = RES11*DJ
RESAB9 = RESAB9*DJ
DRESC = DRESC*DJ
EST = ABS(RES9-RES11)
IF(DRESC.NE.0.0D+00) EST = MAX(EST,DRESC*MIN(1.0D+00,(20.0D+00
* *EST/DRESC)**1.5D+00))
IF(RESAB9.GT.UFLOW) EST = MAX(EMACH*RESAB9,EST)
RETURN
END
subroutine tridv(node,node1,node2,coef,rank)
double precision node(10),node1(10),node2(10),coef
integer rank
double precision s(3),coef1,temp
integer t(3)
coef1=1.0-coef
s(1)=(node(3)-node(5))**2+(node(4)-node(6))**2
s(2)=(node(5)-node(7))**2+(node(6)-node(8))**2
s(3)=(node(3)-node(7))**2+(node(4)-node(8))**2
t(1)=1
t(2)=2
t(3)=3
do 10 i=1,2
do 10 j=i+1,3
if(s(i).lt.s(j)) then
temp=t(i)
t(i)=t(j)
t(j)=temp
end if
10 continue
if(t(rank).eq.1)then
node1(3)=coef*node(3)+coef1*node(5)
node1(4)=coef*node(4)+coef1*node(6)
node1(5)=node(5)
node1(6)=node(6)
node1(7)=node(7)
node1(8)=node(8)
node2(3)=node1(3)
node2(4)=node1(4)
node2(5)=node(7)
node2(6)=node(8)
node2(7)=node(3)
node2(8)=node(4)
else if(t(rank).eq.2) then
node1(3)=coef*node(5)+coef1*node(7)
node1(4)=coef*node(6)+coef1*node(8)
node1(5)=node(7)
node1(6)=node(8)
node1(7)=node(3)
node1(8)=node(4)
node2(3)=node1(3)
node2(4)=node1(4)
node2(5)=node(3)
node2(6)=node(4)
node2(7)=node(5)
node2(8)=node(6)
else
node1(3)=coef*node(3)+coef1*node(7)
node1(4)=coef*node(4)+coef1*node(8)
node1(5)=node(3)
node1(6)=node(4)
node1(7)=node(5)
node1(8)=node(6)
node2(3)=node1(3)
node2(4)=node1(4)
node2(5)=node(5)
node2(6)=node(6)
node2(7)=node(7)
node2(8)=node(8)
end if
node1(9)=coef*node(9)
node2(9)=coef1*node(9)
end
SUBROUTINE INTEGXERROR(MESSG,NMESSG,NERR,LEVEL)
C***BEGIN PROLOGUE XERROR
C***DATE WRITTEN 790801 (YYMMDD)
C***REVISION DATE 820801 (YYMMDD)
C***CATEGORY NO. R3C
C***KEYWORDS ERROR,XERROR PACKAGE
C***AUTHOR JONES, R. E., (SNLA)
C***PURPOSE Processes an error (diagnostic) message.
C***DESCRIPTION
C Abstract
C XERROR processes a diagnostic message, in a manner
C determined by the value of LEVEL and the current value
C of the library error control flag, KONTRL.
C (See subroutine XSETF for details.)
C
C Description of Parameters
C --Input--
C MESSG - the Hollerith message to be processed, containing
C no more than 72 characters.
C NMESSG- the actual number of characters in MESSG.
C NERR - the error number associated with this message.
C NERR must not be zero.
C LEVEL - error category.
C =2 means this is an unconditionally fatal error.
C =1 means this is a recoverable error. (I.e., it is
C non-fatal if XSETF has been appropriately called.)
C =0 means this is a warning message only.
C =-1 means this is a warning message which is to be
C printed at most once, regardless of how many
C times this call is executed.
C
C Examples
C CALL XERROR('SMOOTH -- NUM WAS ZERO.',23,1,2)
C CALL XERROR('INTEG -- LESS THAN FULL ACCURACY ACHIEVED.',
C 43,2,1)
C CALL XERROR('ROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL F
C 1ULLY COLLAPSED.',65,3,0)
C CALL XERROR('EXP -- UNDERFLOWS BEING SET TO ZERO.',39,1,-1)
C
C Latest revision --- 19 MAR 1980
C Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR-
C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES,
C 1982.
C***ROUTINES CALLED XERRWV
C***END PROLOGUE XERROR
CHARACTER*(*) MESSG
C***FIRST EXECUTABLE STATEMENT XERROR
CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.)
RETURN
END
|
