-
Notifications
You must be signed in to change notification settings - Fork 0
/
sepbvp.f90
9350 lines (7826 loc) · 498 KB
/
sepbvp.f90
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
module algprs
use bvpprec
implicit none
public
integer,save :: nminit =8
integer,save :: iprint =0
integer,save :: idum=0
real(dp), save :: uval0 = 0._dp
logical,save :: pdebug =.false.
logical,save :: use_c =.true.
logical,save :: comp_c =.true.
end module algprs
module sepbvp
use bvpprec
implicit none
private
public easybvpc
contains
subroutine easybvpc(nEq,nBc,r,y,alltol)
use algprs, only : uval0
implicit none
integer, intent(in) :: nEq,nBc
real(dp), dimension(:), pointer :: r
real(dp), dimension(:,:), pointer ::y
real(dp), intent(in) :: alltol
!internal working arrays of twpbvpbc
integer, parameter :: nxxdim = 10000000
integer, save :: lwrkfl= 10000000,lwrkin=10000000
integer, parameter :: lwrkMax=100000000
!input of twpbvpbc
logical :: linear, givmsh, giveu
integer :: ncomp,nlbc
integer :: nfxpnt
integer :: ntol
integer :: nmsh,nmax
integer :: nudim, iflbvp
real(dp) :: aleft, aright
real(dp), dimension(1) :: fixpnt
integer, dimension(:), allocatable :: ltol
real(dp), dimension(:), allocatable :: tol
real(dp), dimension(:,:), allocatable :: u
real(dp), dimension(:), allocatable :: xx
real(dp), dimension(:), allocatable :: wrk
integer, dimension(:), allocatable :: iwrk
real(dp) :: ckappa,ckappa1,gamma1
real(dp), dimension(1) :: rpar
integer, dimension(1) :: ipar
!local scope
logical, parameter :: checkInput=.true.
integer :: i,j,n
!interface to the user supplied functions
include 'sepbvp.h'
if ((.not.associated(r)).or.(.not.associated(y))) then
stop 'easybvp: trial solutions not found!'
else
n = size(r,1)
endif
if (checkInput) then
if (n.gt.nxxdim) then
write(*,*)'nmax= ',nxxdim
stop 'r-vector too large'
endif
if ((size(y,1).ne.neq).or.(size(y,2).ne.n)) then
stop 'y-vector size incorrect!'
endif
if ((nBc.gt.neq).or.(nBc.lt.0)) then
write(*,*)'Improper number of boundary conditions!'
write(*,*)'nEq= nBc= ',nEq,nBc
stop
endif
if (r(1).ge.r(n)) stop 'rmin > rmax!'
endif
!boundary points
aleft=r(1)
aright=r(n)
!no fix point added on the mesh
nfxpnt=0
fixpnt(1)=(aleft+aright)/2.
!tolerance vector
ntol=neq
allocate(ltol(ntol),tol(ntol))
do i=1,ntol
ltol(i)=i
tol(i)=alltol
enddo
!mesh and array sizes
nlbc=nBc
ncomp=neq
nudim=ncomp
nmsh=n
!user given mesh
allocate(xx(1:nxxdim))
xx(1:nmsh)=r(1:nmsh)
!the solution vector
allocate(u(nudim,nxxdim))
u(1:ncomp,1:nmsh)=y(1:neq,1:n)
linear = .false.
givmsh=.true.
giveu=.true.
do
allocate(wrk(lwrkfl))
allocate(iwrk(lwrkin))
call twpbvpc(ncomp,nlbc,aleft,aright,nfxpnt,fixpnt,ntol,ltol,tol &
,linear,givmsh,giveu,nmsh,nxxdim,xx,nudim,u,nmax,lwrkfl,wrk &
,lwrkin,iwrk,fsub,dfsub,gsub,dgsub &
,ckappa1,gamma1,ckappa,rpar,ipar,iflbvp)
if (nmax.gt.nxxdim) stop 'nmax > nxxdim'
write(*,*)'twpbvpc on exit: iflbvp= ',iflbvp
write(*,*)'nmsh= ',nmsh
write(*,*)'nmax= ',nmax
select case (iflbvp)
case (0)
if (nmsh.ne.n) then
deallocate(r,y)
allocate(r(nmsh),y(ncomp,nmsh))
r = xx(1:nmsh)
y(1:ncomp,1:nmsh)=u(1:ncomp,1:nmsh)
endif
exit
case (-1)
stop 'easybvp: invalid input parameters!'
case (1)
write(*,*) 'easybvp: maximal numbers of mesh points exceeded!'
!start again with bigger working arrays
write(*,*) '...increasing working array sizes...'
lwrkfl = lwrkfl*5
lwrkin = lwrkin*5
write(*,*) 'lwrkfl= lwrkin= ',lwrkfl,lwrkin
deallocate(wrk,iwrk)
!reset to the trial solution
nlbc=nBc
ncomp=neq
nudim=ncomp
nmsh=n
xx(1:nmsh)=r(1:nmsh)
u(1:ncomp,1:nmsh)=y(1:neq,1:n)
call random_number(uval0)
!check if we are not killing all memory
if (max(lwrkin,lwrkfl).gt.lwrkMax) exit
cycle
case default
stop 'easybvp: internal error!'
end select
enddo
deallocate(ltol,tol,u,xx)
end subroutine easybvpc
!-------------------------------------------------------------------------
!encapsulated into f90 module.
!block data removed + algprs common block converted to algprs module.
!added implicit declaration of integer in between i-o
!external statements replaced by interfaces + include statements.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!credit: http://www.ma.ic.ac.uk/~jcash/BVP_software/readme.php#docs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-------------------------------------------------------------------------
subroutine twpbvpc(ncomp, nlbc, aleft, aright, &
& nfxpnt, fixpnt, ntol, ltol, tol, &
& linear, givmsh, giveu, nmsh, &
& nxxdim, xx, nudim, u, nmax, &
& lwrkfl, wrk, lwrkin, iwrk, &
& fsub, dfsub, gsub, dgsub, &
& ckappa1,gamma1,ckappa,rpar,ipar,iflbvp)
!
! The subroutine twpbvpc is intended to solve two-point boundary
! value problems.
!
!
! References:
! Cash, J. R.; Mazzia, F.
! A new mesh selection algorithm, based on conditioning,
! for two-point boundary value codes.
! J. Comput. Appl. Math. 184 (2005), no. 2, 362--381.
!
! Revision History
!
! revision July 3, 2006
! added rpar and ipar in the functions
! DFSUB, DGSUB, FSUB, GSUB
! changed the name of the variable double in ddouble
!
! revision September 20, 2004
! This is a modified version of twpbvp that uses the conditioning
! in the mesh selection.
!
! New subroutines not included in the old version:
! condestim
! estimkappa
! moncond
! selcondmsh
! selconderrmsh
! smpselcondmsh
!
! Updates subroutines:
! bvpsol
! conv4
! fail4
! fail6
! conv8
! fail8
! newteq
! mshref
! wtchdg
!
! Auxiliary function not used by the old version
! donest
! abdnrm
!
! The common block algprs contains two more variable
! logical pdebug, use_c, comp_c
! common/algprs/ nminit, pdebug, iprint, idum, uval0, use_c, comp_c
!
! The starting subroutine is twpbvp by J.R. Cash and M.H. Wright
!
!module
use algprs
implicit integer (i-n)
implicit real(dp) (a-h,o-z)
dimension rpar(*),ipar(*)
dimension fixpnt(*), ltol(*), tol(*)
dimension xx(*), u(nudim,*)
dimension wrk(lwrkfl), iwrk(lwrkin)
logical linear, givmsh, giveu
! external fsub, dfsub, gsub, dgsub
! logical pdebug, use_c, comp_c
! common/algprs/ nminit, pdebug, iprint, idum, uval0, use_c, comp_c
intrinsic abs, min
parameter ( zero = 0.0_dp )
!interface
include 'sepbvp.h'
! Check for invalid input parameters. If any parameters are
! invalid, exit with the flag iflbvp set to -1.
iflbvp = -1
if (ncomp .le. 0) return
if (nlbc .lt. 0 .or. nlbc .gt. ncomp) return
if (aleft .ge. aright) return
if (nfxpnt .lt. 0) return
if (givmsh .and. nmsh .lt. nfxpnt+2) return
if (givmsh .and. xx(1) .ne. aleft) return
! SCMODIFIED add an extra condition to avoid accessing xx(0)
if (nmsh .gt. 0) then
if (givmsh .and. xx(nmsh) .ne. aright) return
end if
if (nfxpnt .gt. 0) then
if (fixpnt(1) .le. aleft) return
if (fixpnt(nfxpnt) .ge. aright) return
do 50 i = 1, nfxpnt-1
if (fixpnt(i+1) .le. fixpnt(i)) return
50 continue
endif
if (ntol .lt. 1) return
do 60 i = 1, ntol
if (ltol(i) .lt. 0 .or. ltol(i) .gt. ncomp) return
if (tol(i) .le. zero) return
60 continue
if (giveu .and. .not. givmsh) return
if (use_c .and. .not. comp_c) return
if (nudim .le. 0) return
if (lwrkfl .le. 0 .or. lwrkin .le. 0) return
! Calculate maximum number of mesh points possible with the
! given floating-point and integer workspace.
isp = lwrkfl - 2*ntol - 22*ncomp - 6*ncomp*ncomp
iden = 6*ncomp*ncomp + 22*ncomp + 3
nmax1 = isp/iden
isp = lwrkin - 2*ncomp
nmax2 = isp/(2*ncomp+3)
nmax = min(nmax1,nmax2)
! nmax from workspace
nmax = min(nmax, nxxdim)
! nmax from size of u and xx
if (iprint .ge. 0) write(6,901) nmax
901 format(1h ,'nmax from workspace =',i8)
if (nmax .le. 1) return
! Partition floating point workspace.
irhs = 1
lrhs = ncomp*nmax
itpblk = irhs + lrhs
ltpblk = ncomp*nlbc
ibtblk = itpblk + ltpblk
lbtblk = ncomp*(ncomp - Nlbc)
iajac = ibtblk + lbtblk
lajac = 2*ncomp*ncomp*nmax
ibhold = iajac + lajac
lbhold = ncomp*ncomp*nmax
ichold = ibhold + lbhold
lchold = ncomp*ncomp*nmax
ifval = ichold + lchold
lfval = ncomp*nmax
idef = ifval + lfval
ldef = ncomp*(nmax-1)
idefex = idef + ldef
ldefex = ncomp*(nmax-1)
! def6 uses the same space as defexp
idef6 = idefex
ldef6 = ncomp*(nmax-1)
idefim = idef6 + ldef6
ldefim = ncomp*(nmax-1)
! def8 uses the same space as defimp
idef8 = idefim
ldef8 = ncomp*(nmax-1)
iusve = idef8 + ldef8
lusve = ncomp*nmax
iuold = iusve + lusve
luold = ncomp*nmax
itmrhs = iuold + luold
ltmrhs = ncomp*nmax
irhtri = itmrhs + ltmrhs
lrhtri = ncomp*nmax
idelu = irhtri + lrhtri
ldelu = ncomp*nmax
ixmer = idelu + ldelu
lxmer = ncomp*nmax
! rerr occupies the same space as xmerit
irerr = ixmer
lrerr = ncomp*nmax
iutri = irerr + lrerr
lutri = ncomp*nmax
iermx = iutri + lutri
lermx = nmax
irtdc = iermx + lermx
lrtdc = nmax
ixxold = irtdc + lrtdc
lxxold = nmax
iuint = ixxold + lxxold
luint = ncomp
iftmp = iuint + luint
lftmp = ncomp
idgtm = iftmp + lftmp
ldgtm = ncomp
idftm1 = idgtm + ldgtm
ldftm1 = ncomp*ncomp
idftm2 = idftm1 + ldftm1
ldftm2 = ncomp*ncomp
itmp = idftm2 + ldftm2
ltmp = ncomp*8
idsq = itmp + ltmp
ldsq = ncomp*ncomp
idexr = idsq + ldsq
ldexr = ncomp
ietst6 = idexr + ldexr
letst6 = ntol
ietst8 = ietst6 + letst6
letst8 = ntol
iamg = ietst8 + letst8
lamg = ncomp*nmax
ic1 = iamg + lamg
lc1 = ncomp*ncomp*nmax
iwrkrhs = ic1+lc1
lwrkrhs = ncomp*nmax
ilast = iwrkrhs + lwrkrhs
if (iprint .eq. 1) write(6,903) ilast
903 format(1h ,'ilast',i10)
! Partition integer workspace.
iiref = 1
liref = nmax
iihcom = iiref + liref
lihcom = nmax
iipvbk = iihcom + lihcom
lipvbk = ncomp*nmax
iipvlu = iipvbk + lipvbk
lipvlu = ncomp
iisign = iipvlu + lipvlu
lisign = ncomp*nmax
call bvpsol(ncomp, nmsh, nlbc, aleft, aright, &
& nfxpnt, fixpnt, ntol, ltol, tol, nmax, linear, &
& giveu, givmsh, xx, nudim, u, &
& wrk(idefex), wrk(idefim), wrk(idef), wrk(idelu), &
& wrk(irhs), wrk(ifval), &
& wrk(itpblk), wrk(ibtblk), wrk(iajac), wrk(ibhold), &
& wrk(ichold), iwrk(iipvbk), iwrk(iipvlu),iwrk(iisign), &
& wrk(iuint), wrk(iftmp), wrk(itmrhs), &
& wrk(idftm1), wrk(idftm2), wrk(idgtm), &
& wrk(iutri), wrk(irhtri), wrk(ixmer), &
& wrk(ixxold), wrk(iuold), wrk(iusve), &
& wrk(itmp), wrk(idsq), wrk(idexr), wrk(irtdc), &
& wrk(irerr), wrk(ietst6), wrk(ietst8), wrk(iermx), &
& iwrk(iihcom), iwrk(iiref), wrk(idef6), wrk(idef8), &
& fsub,dfsub,gsub,dgsub,iflbvp, &
& wrk(iamg),wrk(ic1),wrk(iwrkrhs), &
& ckappa1,gamma1,ckappa,rpar, ipar)
return
end subroutine
subroutine bvpsol(ncomp, nmsh, nlbc, aleft, aright, &
& nfxpnt, fixpnt, &
& ntol, ltol, tol, nmax, linear, giveu, givmsh, &
& xx, nudim, u, defexp, defimp, def, delu, rhs, fval, &
& topblk, botblk, ajac, bhold, chold, ipvblk, ipivlu,isign, &
& uint, ftmp, tmprhs, dftmp1, dftmp2, dgtm, &
& utrial, rhstri, xmerit, xxold, uold, usave, &
& tmp, dsq, dexr, ratdc, rerr, &
& etest6, etest8, ermx, ihcomp, irefin, &
& def6, def8, fsub, dfsub, gsub, dgsub, iflbvp, &
& amg, c1, wrkrhs,ckappa1,gamma1,ckappa,rpar,ipar)
!module
use algprs
implicit integer (i-n)
implicit real(dp) (a-h,o-z)
dimension rpar(*), ipar(*)
dimension fixpnt(*), ltol(ntol), tol(ntol)
dimension xx(*), u(nudim, *)
dimension defexp(ncomp,*), defimp(ncomp,*), def(ncomp,*)
dimension delu(ncomp, *), rhs(*), fval(ncomp,*)
dimension topblk(nlbc,*), botblk(ncomp-nlbc,*)
dimension ajac(ncomp, 2*ncomp, *)
dimension bhold(ncomp, ncomp, *), chold(ncomp, ncomp, *)
dimension ipivlu(*), ipvblk(*), isign(*)
dimension uint(ncomp), ftmp(ncomp)
dimension dgtm(ncomp), tmprhs(*)
dimension dftmp1(ncomp, ncomp), dftmp2(ncomp, ncomp)
dimension utrial(ncomp,*), rhstri(*)
dimension xmerit(ncomp, *)
dimension xxold(*), uold(ncomp,*), usave(ncomp,*)
dimension tmp(ncomp,8)
dimension dsq(ncomp,ncomp), dexr(ncomp)
dimension ratdc(*), rerr(ncomp,*)
dimension etest6(*), etest8(*), ermx(*)
dimension ihcomp(*), irefin(*)
dimension def6(ncomp,*), def8(ncomp,*)
dimension amg(*), c1(ncomp,*), wrkrhs(*)
logical linear, giveu, givmsh, ddouble
! external fsub, dfsub, gsub, dgsub
common/mchprs/flmin, flmax, epsmch
! logical pdebug, use_c, comp_c
! common/algprs/ nminit, pdebug, iprint, idum, uval0, use_c, comp_c
intrinsic max
logical smooth, succes, strctr, trst6, reaft6
logical onto6, onto8, ludone, rhsgiv, maxmsh
logical first4, first8
logical mchset
save mchset
logical frscal
! save frscal
logical stab_kappa, stab_gamma, stab_cond, stiff_cond
parameter (zero = 0.0_dp, one = 1.0_dp)
parameter (third = 0.33_dp, fourth = 0.25_dp)
parameter (quan6 = 0.1_dp )
integer, parameter ::itcondmax = 10
! blas: dload
! real(dp) d1mach
!interface
include 'sepbvp.h'
data mchset/.true./
data fxfct/10.0_dp/
data maxmsh/.false./
frscal = .true.
if (mchset) then
flmin = d1mach(1)
flmax = d1mach(2)
epsmch = d1mach(3)
if (pdebug) write(6,901) epsmch
mchset = .false.
endif
! The routine stcons calculates integration constants stored in
! labeled common consts.
call stcons
! Set up arrays for the error tests.
if (.not. linear) then
call dload(ntol, one, etest6, 1)
else
do 10 i = 1, ntol
etest6(i) = one/max(quan6, tol(i)**third)
10 continue
endif
nmold = 1
smooth = .false.
strctr = .false.
trst6 = .true.
reaft6 = .false.
numbig = 0
nummed = 0
first4 = .true.
first8 = .true.
onto6 = .false.
! initialize parameter for the conditioning estimation
itcond = 0
gamma1old = flmax
gamma1 = flmax
ckappa1old = flmax
ckappa1 = flmax
ckappa = flmax
stiff_cond = .false.
stab_cond = .false.
tolmin = flmax
do i=1,ntol
tolmin = min(tol(i),tolmin)
end do
! If givmsh is .true., the initial number of mesh points must be
! provided by the user in nmsh, and the mesh points must be
! contained in the array xx (of dimension nmsh).
! Otherwise, nmsh is set to its default value, and a
! uniform initial mesh is created.
if (.not. giveu .and. .not. givmsh) then
nmsh = nminit
if (nmsh .lt. nfxpnt+2) nmsh = nfxpnt + 2
call unimsh(nmsh, aleft, aright, nfxpnt, fixpnt, xx)
endif
if (pdebug) then
write(6,902)
call sprt(nmsh, xx)
endif
if (.not. giveu) call initu(ncomp, nmsh, xx, nudim, u, &
& rpar,ipar )
!**** top of logic for 4th order solution ****
400 continue
if (iprint .eq. 1) write(6,903) nmsh
! Set the def (deferred correction) array to zero.
call mtload(ncomp, nmsh-1, zero, ncomp, def)
iorder = 4
! The routine fneval calls fsub at the mesh xx and the
! solution u, and saves the values in the array fval.
call fneval(ncomp, nmsh, xx, nudim, u, fval, fsub,rpar,ipar)
! Try to compute a 4th order solution by solving a system of nonlinear
! equations.
if (linear) then
ludone = .false.
call lineq( ncomp, nmsh, nlbc, &
& ludone, xx, nudim, u, def, &
& delu, rhs, fval, &
& uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, &
& ajac, topblk, botblk, bhold, chold, ipvblk, &
& fsub, dfsub, gsub, dgsub, iflnwt,rpar,ipar)
! Call fneval to evaluate the fval array at the new solution u.
! (Such a call is not necessary for the nonlinear case because
! fval is called within newteq for the new u.)
call fneval(ncomp, nmsh, xx, nudim, u, fval, fsub,rpar,ipar)
else
rhsgiv = .false.
call newteq(ncomp, nmsh, nlbc, &
& rhsgiv, ntol, ltol, tol, &
& xx, nudim, u, def, &
& delu, rhs, fval, &
& utrial, rhstri, &
& uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, xmerit, &
& ajac, topblk, botblk, bhold, chold, ipvblk, &
& fsub, dfsub, gsub, dgsub, itnwt, iflnwt,rpar,ipar,frscal)
endif
!
! these flags are used in the mesh selection strategy
!
if (iflnwt .eq. 0) then
! COMPUTE ESTIMATIONS OF CONDITION NUMBERS KAPPA1 and GAMMA1
! related to perturbation to the boundary conditions
!
N =nmsh*ncomp
ninter=nmsh-1
if (comp_c) then
gamma1old = gamma1
ckappa1old = ckappa1
call CONDESTIM(aleft,aright,nmsh,ncomp,N,xx,topblk,nlbc, &
& ncomp, ajac, ncomp,2*ncomp,ninter,botblk,ncomp-nlbc, &
& ipvblk,isign,amg,c1,wrkrhs,ckappa1,gamma1)
if (iprint .ge. 0) then
!
! COMPUTE ESTIMATION OF THE CONDITION NUMBER KAPPA
!
call ESTIMKAPPA(nmsh,ncomp,N,xx,topblk, &
& nlbc,ncomp,ajac, ncomp,2*ncomp, &
& ninter,botblk,ncomp-nlbc,ipvblk,isign,c1,wrkrhs,ckappa)
end if
if (iprint .ge. 0) then
write(6,1001) ckappa1/gamma1
write(6,1002) gamma1
write(6,1003) ckappa1
write(6,1004) ckappa
end if
stab_kappa = abs(ckappa1old-ckappa1)/(1d0+ckappa1).lt.5d-2 &
& .and. ckappa1 .lt. flmax .and. gamma1.lt.flmax
stab_gamma = abs(gamma1old-gamma1)/(1d0+gamma1).lt.5d-2 &
& .and. gamma1.lt.flmax .and. ckappa1 .lt. flmax
stab_cond = stab_kappa .and. stab_gamma
stiff_cond = (( (ckappa1/gamma1 .ge. 10d0 )))
if (iprint .eq. 1) then
write(6,*) 'stab_kappa = ', stab_kappa
write(6,*) 'stab_gamma = ', stab_gamma
write(6,*) 'stiff_cond = ', stiff_cond
end if
end if
! endif if (comp_c)
call conv4( ncomp, nmsh, ntol, ltol, tol, linear, nmax, &
& xx, nudim, u, defexp, defimp, def, fval, &
& tmp, bhold, chold, dsq, dexr, usave, &
& ratdc, rerr, ipivlu, nmold, xxold, &
& smooth, reaft6, onto6, strctr, trst6, ddouble , &
& fsub, maxmsh, succes, first4, &
& amg,stab_cond,ckappa,stiff_cond,wrkrhs, &
& nfxpnt, fixpnt, irefin,rpar,ipar)
if (pdebug .and. .not. onto6) write (6,904)
else
if (comp_c) then
if (iflnwt .ne. -1) then
gamma1old = gamma1
ckappa1old = ckappa1
N =nmsh*ncomp
ninter=nmsh-1
call CONDESTIM(aleft,aright,nmsh,ncomp,N,xx,topblk,nlbc, &
& ncomp, ajac, ncomp,2*ncomp,ninter,botblk,ncomp-nlbc, &
& ipvblk,isign,amg,c1,wrkrhs,ckappa1,gamma1)
if (iprint .ge.0) then
call ESTIMKAPPA(nmsh,ncomp,N,xx,topblk, &
& nlbc,ncomp,ajac, ncomp,2*ncomp, &
& ninter,botblk,ncomp-nlbc,ipvblk,isign,c1,wrkrhs,ckappa)
end if
end if
if (iprint .ge. 0) then
write(6,1001) ckappa1/gamma1
write(6,1002) gamma1
write(6,1003) ckappa1
write(6,1004) ckappa
end if
stab_kappa = abs(ckappa1old-ckappa1)/(1d0+ckappa1).lt.5d-2 &
& .and. ckappa1 .lt. flmax .and. gamma1.lt.flmax
stab_gamma = abs(gamma1old-gamma1)/(1d0+gamma1).lt.5d-2 &
& .and. gamma1.lt.flmax .and. ckappa1 .lt. flmax
stab_cond = stab_kappa .and. stab_gamma
stiff_cond = (( (ckappa1/gamma1 .ge. 5d0 )))
if (iprint .eq. 1) then
write(6,*) 'stab_kappa = ',stab_kappa
write(6,*) 'stab_gamma = ', stab_gamma
write(6,*) 'stiff_cond = ', stiff_cond
end if
end if
! end if if(comp_c)
succes = .false.
onto6 = .false.
reaft6 = .false.
call fail4( ncomp, nmsh, nlbc, ntol, ltol, &
& xx, nudim, u, rhs, linear, nmax, &
& nmold, xxold, uold, ratdc, &
& iorder, iflnwt, itnwt, ddouble , maxmsh, &
& numbig, nummed,wrkrhs,amg,stab_cond,stiff_cond, &
& nfxpnt, fixpnt, irefin,itcond,itcondmax,rpar,ipar)
endif
if (succes) then
if (comp_c) then
call ESTIMKAPPA(nmsh,ncomp,N,xx,topblk, &
& nlbc,ncomp,ajac, ncomp,2*ncomp, &
& ninter,botblk,ncomp-nlbc,ipvblk,isign,c1,wrkrhs,ckappa)
end if
if (iprint .ne. -1 .and. comp_c ) then
if (ckappa .ge. tolmin/epsmch) write(6,1005)
end if
iflbvp = 0
return
elseif (maxmsh) then
go to 900
elseif (.not. onto6) then
go to 400
endif
! To reach here, onto6 must be .true.
!*** logic for 6th order ****
if (iprint .eq. 1) write(6,905)
! Save the 4th order solution on this mesh in uold.
call matcop(nudim, ncomp, ncomp, nmsh, u, uold)
! Copy the current mesh into the xxold array.
nmold = nmsh
call dcopy(nmold, xx, 1, xxold, 1)
iorder = 6
if (linear) then
call lineq( ncomp, nmsh, nlbc, &
& ludone, xx, nudim, u, def, &
& delu, rhs, fval, &
& uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, &
& ajac, topblk, botblk, bhold, chold, ipvblk, &
& fsub, dfsub, gsub, dgsub, iflnwt,rpar,ipar)
else
call fixjac(ncomp, nmsh, nlbc, &
& iorder, ntol, ltol, tol, &
& xx, nudim, u, def, def, delu, &
& rhs, fval, utrial, rhstri, &
& rnsq, uint, ftmp, tmprhs, &
& ajac, topblk, botblk, ipvblk, &
& fsub, gsub, iflnwt,rpar,ipar)
! If the fixed Jacobian iterations fail but rnsq is small,
! try a Newton procedure. Set rhsgiv to indicate that
! the right-hand side and fval have already been evaluated
! at the given u.
if (iflnwt .eq. -3 .and. rnsq .lt. fxfct*epsmch) then
rhsgiv = .true.
call newteq(ncomp, nmsh, nlbc, &
& rhsgiv, ntol, ltol, tol, &
& xx, nudim, u, def, &
& delu, rhs, fval, &
& utrial, rhstri, &
& uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, xmerit, &
& ajac, topblk, botblk, bhold, chold, ipvblk, &
& fsub, dfsub, gsub, dgsub, iter, iflnwt,rpar,ipar,frscal)
endif
endif
if (iflnwt .eq. 0) then
call conv6(ncomp, nmsh, ntol, ltol, tol, &
& nudim, u, uold, etest6, err6, &
& trst6, onto8, reaft6, succes)
else
onto8 = .false.
call fail6( ncomp, nmsh, nlbc, ntol, ltol, tol, &
& nfxpnt, fixpnt, iorder, nmax, &
& xx, nudim, u, rhs, usave, xxold, uold, nmold, &
& ihcomp, irefin, &
& rerr, ermx, ratdc, &
& reaft6, ddouble , succes, maxmsh, &
& numbig, nummed, &
& wrkrhs,amg, stab_cond,ckappa1,gamma1,ckappa, &
& stiff_cond, itcond, itcondmax)
endif
if (succes) then
if (comp_c) then
call ESTIMKAPPA(nmsh,ncomp,N,xx,topblk, &
& nlbc,ncomp,ajac, ncomp,2*ncomp, &
& ninter,botblk,ncomp-nlbc,ipvblk,isign,c1,wrkrhs,ckappa)
end if
if (iprint .ne. -1 .and. comp_c ) then
if ( ckappa .ge. tolmin/epsmch) write(6,1005)
end if
iflbvp = 0
return
elseif (maxmsh) then
go to 900
elseif (.not. onto8) then
go to 400
endif
!**** logic for trying to calculate 8th order solution *****
if (iprint .eq. 1) write(6,906)
call matcop(nudim, ncomp, ncomp, nmsh, u, uold)
! Copy the current mesh into the xxold array.
nmold = nmsh
call dcopy(nmold, xx, 1, xxold, 1)
! Save the old deferred correction vector def in def6.
call matcop(ncomp, ncomp, ncomp, nmsh-1, def, def6)
! For linear problems, calculate the fval array for the
! new solution u.
if (linear) call fneval(ncomp, nmsh, xx, nudim, u, fval, fsub, &
& rpar,ipar)
! Calculate 8th order deferred corrections (the array def8).
call df8cal (ncomp, nmsh, xx, nudim, u, fval, def8, &
& tmp, fsub,rpar,ipar)
! For linear problems, the def array is the def8 array.
! For nonlinear problems, add the def8 array to the
! already-calculated def array.
if (linear) then
call matcop(ncomp, ncomp, ncomp, nmsh-1, def8, def)
else
call maxpy(ncomp, nmsh-1, one, def8, ncomp, def)
endif
iorder = 8
if (linear) then
call lineq( ncomp, nmsh, nlbc, &
& ludone, xx, nudim, u, def, &
& delu, rhs, fval, &
& uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, &
& ajac, topblk, botblk, bhold, chold, ipvblk, &
& fsub, dfsub, gsub, dgsub, iflnwt,rpar,ipar)
else
call fixjac(ncomp, nmsh, nlbc, &
& iorder, ntol, ltol, tol, &
& xx, nudim, u, def, def8, delu, &
& rhs, fval, utrial, rhstri, &
& rnsq, uint, ftmp, tmprhs, &
& ajac, topblk, botblk, ipvblk, &
& fsub, gsub, iflnwt,rpar,ipar)
! If the fixed Jacobian iterations fail but rnsq is small,
! try a Newton procedure. Set rhsgiv to indicate that
! the right-hand side and fval have already been evaluated