forked from geodynamics/seismic_cpml
-
Notifications
You must be signed in to change notification settings - Fork 0
/
attenuation_model_with_SolvOpt.f90
2000 lines (1801 loc) · 59.4 KB
/
attenuation_model_with_SolvOpt.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
! use of SolvOpt to compute attenuation relaxation mechanisms,
! from Emilie Blanc, Bruno Lombard and Dimitri Komatitsch, CNRS Marseille, France, for a Generalized Zener body model.
! The SolvOpt algorithm was developed by Franz Kappel and Alexei V. Kuntsevich
! and is available open source at http://www.uni-graz.at/imawww/kuntsevich/solvopt
!
! It is described in Kappel and Kuntsevich, An Implementation of Shor's r-Algorithm,
! Computational Optimization and Applications, vol. 15, p. 193-205 (2000).
! If you use this code for your own research, please cite some (or all) of these articles:
!
! @Article{BlKoChLoXi15,
! Title = {Positivity-preserving highly-accurate optimization of the {Z}ener viscoelastic model, with application
! to wave propagation in the presence of strong attenuation},
! Author = {\'Emilie Blanc and Dimitri Komatitsch and Emmanuel Chaljub and Bruno Lombard and Zhinan Xie},
! Journal = {Geophysical Journal International},
! Year = {2015},
! Note = {in press.}}
!-------------------------------------------------------------------------
! From Bruno Lombard, May 2014:
! En interne dans le code ci-dessous on travaille en (Theta, Kappa).
! Les Theta sont les points et les Kappa sont les poids.
! Pour repasser en (Tau_Sigma, Tau_Epsilon), on doit appliquer les formules:
!
! Tau_Sigma = 1 / Theta
! Tau_Epsilon = (1 / Theta) * (1 + Nrelax * Kappa) = Tau_Sigma * (1 + Nrelax * Kappa)
! The system to solve can be found in equation (7) of:
! Lombard and Piraux, Numerical modeling of transient two-dimensional viscoelastic waves,
! Journal of Computational Physics, Volume 230, Issue 15, Pages 6099-6114 (2011)
! Suivant les compilateurs et les options de compilation utilisees,
! il peut y avoir des differences au 4eme chiffre significatif. C'est sans consequences sur la precision du calcul :
! l'erreur est de 0.015 % avec optimisation non lineaire, a comparer a 1.47 % avec Emmerich and Korn (1987).
! Si je relance le calcul en initialisant avec le resultat precedent, ce chiffre varie a nouveau tres legerement.
!-------------------------------------------------------------------------
! From Bruno Lombard, June 2014:
! j'ai relu en detail :
! [1] Carcione, Kosslof, Kosslof, "Viscoacoustic wave propagation simulation in the Earth",
! Geophysics 53-6 (1988), 769-777
!
! [2] Carcione, Kosslof, Kosslof, "Wave propagation simulation in a linear viscoelastic medium",
! Geophysical Journal International 95 (1988), 597-611
!
! [3] Moczo, Kristek, "On the rheological models used for time-domain methods of seismic wave propagation",
! Geophysical Research Letters 32 (2005).
! Le probleme provient probablement d'une erreur recurrente dans [1,2] et datant de Liu et al 1976 :
! l'oubli du facteur 1/N dans la fonction de relaxation d'un modele de Zener a N elements.
! Il est effectivement facile de faire l'erreur. Voir l'equation (12) de [3], et le paragraphe qui suit.
! Du coup le module de viscoelasticite est faux dans [1,2], et donc le facteur de qualite,
! et donc les temps de relaxation tau_sigma...
! Apres, [2] calcule une solution analytique juste, mais avec des coefficients sans sens physique.
! Et quand SPECFEM2D obtient un bon accord avec cette solution analytique, ca valide SPECFEM, mais pas le calcul des coefficients.
! Il y a donc une erreur dans [1,2], et [3] a raison.
! Sa solution analytique decoule d'un travail sur ses fonctions de relaxation (A4),
! qu'il injecte ensuite dans la relation de comportement (A1) et travaille en Fourier.
! Le probleme est que sa fonction de relaxation (A4) est fausse : il manque 1/N.
! De ce fait, sa solution analytique est coherente avec sa solution numerique.
! Dans les deux cas, ce sont les memes temps de relaxation qui sont utilises. Mais ces temps sont calcules de facon erronee.
!-------------------------------------------------------------------------
! From Dimitri Komatitsch, June 2014:
! In [2] Carcione, Kosslof, Kosslof, "Wave propagation simulation in a linear viscoelastic medium",
! Geophysical Journal International 95 (1988), 597-611
! there is another mistake: in Appendix B page 611 Carcione writes omega/(r*v),
! but that is not correct, it should be omega*r/v instead.
!---------------------------------------------------
! From Emilie Blanc, April 2014:
! le programme SolvOpt d'optimisation non-lineaire
! avec contrainte. Ce programme prend quatre fonctions en entree :
! - fun() est la fonction a minimiser
! - grad() est le gradient de la fonction a minimiser par rapport a chaque parametre
! - func() est le maximum des residus (= 0 si toutes les contraintes sont satisfaites)
! - gradc() est le gradient du maximum des residus (= 0 si toutes les
! contraintes sont satisfaites)
! Ce programme a ete developpe par Kappel et Kuntsevich. Leur article decrit l'algorithme.
! J'ai utilise ce code pour la poroelasticite haute-frequence, et aussi en
! viscoelasticite fractionnaire (modele d'Andrade, avec Bruno Lombard et
! Cedric Bellis). Nous pouvons interagir sur l'algorithme d'optimisation
! pour votre modele visco, et etudier l'effet des coefficients ainsi obtenus.
!---------------------------------------------------
! From Emilie Blanc, March 2014:
! Les entrees du programme principal sont le nombre de variables
! diffusives, le facteur de qualite voulu Qref et la frequence centrale f0.
! Cependant, pour l'optimisation non-lineaire, j'ai mis theta_max=100*f0
! et non pas theta_max=2*pi*100*f0. En effet, dans le programme, on
! travaille sur les frequences, et non pas sur les frequences angulaires.
! Cela dit, dans les deux cas j'obtiens les memes coefficients...
!---------------------------------------------------
SUBROUTINE compute_attenuation_coeffs(N,Qref,f0,f_min,f_max,tau_epsilon,tau_sigma)
IMPLICIT NONE
! pi
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: TWO_PI = 2.d0 * PI
integer, intent(in) :: N
double precision, intent(in) :: Qref,f_min,f_max,f0
double precision, dimension(1:N), intent(out) :: tau_epsilon,tau_sigma
integer i
double precision, dimension(1:N) :: point,weight
! nonlinear optimization with constraints
call nonlinear_optimization(N,Qref,f0,point,weight,f_min,f_max)
do i = 1,N
tau_sigma(i) = 1.d0 / point(i)
tau_epsilon(i) = tau_sigma(i) * (1.d0 + N * weight(i))
enddo
! print *,'points = '
! do i = 1,N
! print *,point(i)
! enddo
! print *
! print *,'weights = '
! do i = 1,N
! print *,weight(i)
! enddo
! print *
! print *,'tau_epsilon = '
! do i = 1,N
! print *,tau_epsilon(i)
! enddo
! print *
! print *,'tau_sigma = '
! do i = 1,N
! print *,tau_sigma(i)
! enddo
! print *
END SUBROUTINE compute_attenuation_coeffs
!---------------------------------------------------
! classical calculation of the coefficients based on linear least squares
SUBROUTINE decomposition_LU(a,i_min,n,indx,d)
IMPLICIT NONE
integer, intent(in) :: i_min,n
double precision, intent(out) :: d
integer, dimension(i_min:n), intent(inout) :: indx
double precision, dimension(i_min:n,i_min:n), intent(inout) :: a
integer i,imax,j,k
double precision big,dum,somme,eps
double precision, dimension(i_min:n) :: vv
imax = 0
d = 1.
eps = 1.e-20
do i = i_min,n
big = 0.
do j = i_min,n
if (abs(a(i,j))>big) then
big = abs(a(i,j))
endif
enddo
if (big==0.) then
print*,'Singular matrix in routine decomposition_LU'
endif
vv(i) = 1./big
enddo
do j = i_min,n
do i = i_min,j-1
somme = a(i,j)
do k = i_min,i-1
somme = somme - a(i,k)*a(k,j)
enddo
a(i,j) = somme
enddo
big = 0.
do i = j,n
somme = a(i,j)
do k = i_min,j-1
somme = somme - a(i,k)*a(k,j)
enddo
a(i,j) = somme
dum = vv(i)*abs(somme)
if (dum>=big) then
big = dum
imax = i
endif
enddo
if (j/=imax) then
do k = i_min,n
dum = a(imax,k)
a(imax,k) = a(j,k)
a(j,k) = dum
enddo
d = -d
vv(imax) = vv(j)
endif
indx(j) = imax
if (a(j,j) == 0.) then
a(j,j) = eps
endif
if (j/=n) then
dum = 1./a(j,j)
do i = j+1,n
a(i,j) = a(i,j)*dum
enddo
endif
enddo
END SUBROUTINE decomposition_LU
SUBROUTINE LUbksb(a,i_min,n,indx,b,m)
IMPLICIT NONE
integer, intent(in) :: i_min,n,m
integer, dimension(i_min:n), intent(in) :: indx
double precision, dimension(i_min:n,i_min:n), intent(in) :: a
double precision, dimension(i_min:n,i_min:m), intent(inout) :: b
integer i,ip,j,ii,k
double precision somme
do k = i_min,m
ii = -1
do i = i_min,n
ip = indx(i)
somme = b(ip,k)
b(ip,k) = b(i,k)
if (ii/=-1) then
do j = ii,i-1
somme = somme - a(i,j)*b(j,k)
enddo
else if (somme/=0.) then
ii = i
endif
b(i,k) = somme
enddo
do i = n,i_min,-1
somme = b(i,k)
do j = i+1,n
somme = somme - a(i,j)*b(j,k)
enddo
b(i,k) = somme/a(i,i)
enddo
enddo
END SUBROUTINE LUbksb
SUBROUTINE syst_LU(a,i_min,n,b,m)
IMPLICIT NONE
integer, intent(in) :: i_min,n,m
double precision, dimension(i_min:n,i_min:n), intent(in) :: a
double precision, dimension(i_min:n,i_min:m), intent(inout) :: b
integer i,j
integer, dimension(i_min:n) :: indx
double precision d
double precision, dimension(i_min:n,i_min:n) :: aux
do j = i_min,n
indx(j) = 0
do i = i_min,n
aux(i,j) = a(i,j)
enddo
enddo
call decomposition_LU(aux,i_min,n,indx,d)
call LUbksb(aux,i_min,n,indx,b,m)
END SUBROUTINE syst_LU
SUBROUTINE lfit_zener(x,y,sig,ndat,poids,ia,covar,chisq,ma,Qref,point)
! ma = nombre de variable diffusive
! ndat = m = K nombre d'abcisse freq_k
IMPLICIT NONE
integer, intent(in) :: ndat,ma
logical, dimension(1:ma), intent(in) :: ia
double precision, intent(in) :: Qref
double precision, intent(out) :: chisq
double precision, dimension(1:ndat), intent(in) :: x,y,sig
double precision, dimension(1:ma), intent(in) :: point
double precision, dimension(1:ma), intent(out) :: poids
double precision, dimension(1:ma,1:ma), intent(out) :: covar
integer i,j,k,l,mfit
double precision ym,wt,sig2i
double precision, dimension(1:ma) :: afunc
double precision, dimension(1:ma,1:1) :: beta
mfit = 0
do j = 1,ma
if(ia(j)) then
mfit = mfit + 1
endif
enddo
if(mfit==0) then
print*,'lfit: no parameters to be fitted'
endif
do j=1,mfit
beta(j,1) = 0.
do k=1,mfit
covar(j,k) = 0.
enddo
enddo
do i=1,ndat
call func_zener(x(i),afunc,ma,Qref,point)
ym = y(i)
if (mfit < ma) then
do j=1,ma
if (.NOT. ia(j)) then
ym = ym - poids(j) * afunc(j)
endif
enddo
endif
sig2i = 1. / (sig(i) * sig(i))
j = 0
do l=1,ma
if (ia(l)) then
j = j+1
wt = afunc(l) * sig2i
k = count(ia(1:l))
covar(j,1:k) = covar(j,1:k) + wt * pack(afunc(1:l),ia(1:l))
beta(j,1) = beta(j,1) + ym * wt
endif
enddo
enddo
do j=2,mfit,1
do k=1,j-1,1
covar(k,j) = covar(j,k)
enddo
enddo
if (ma==1) then
poids(1) = beta(1,1)/covar(1,1)
else if (ma>1) then
call syst_LU(covar,1,mfit,beta,1)
poids(1:ma) = unpack(beta(1:ma,1),ia,poids(1:ma))
endif
chisq = 0.
do i=1,ndat
call func_zener(x(i),afunc,ma,Qref,point)
chisq=chisq+((y(i)-dot_product(poids(1:ma),afunc(1:ma)))/sig(i))**2
enddo
END SUBROUTINE lfit_zener
SUBROUTINE func_zener(x,afunc,N,Qref,point)
IMPLICIT NONE
integer, intent(in) :: N
double precision, intent(in) :: x,Qref
double precision, dimension(1:N), intent(in) :: point
double precision, dimension(1:N), intent(out) :: afunc
integer k
double precision num,deno
do k = 1,N
num = x * (point(k) - x / Qref)
deno = point(k) * point(k) + x * x
afunc(k) = num / deno
enddo
END SUBROUTINE func_zener
SUBROUTINE remplit_point(fmin,fmax,N,point)
IMPLICIT NONE
! pi
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: TWO_PI = 2.d0 * PI
integer, intent(in) :: N
double precision, intent(in) :: fmin,fmax
double precision, dimension(1:N), intent(out) :: point
integer l
IF (N == 1) THEN
point(1) = sqrt(fmin * fmax)
ELSE
do l = 1, N, 1
point(l) = (fmax/fmin) ** ((l-1.)/(N-1.))
point(l) = TWO_PI * point(l) * fmin
enddo
endif
END SUBROUTINE remplit_point
SUBROUTINE classical_linear_least_squares(Qref,poids,point,N,fmin,fmax)
IMPLICIT NONE
! pi
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: TWO_PI = 2.d0 * PI
integer, intent(in) :: N
double precision, intent(in) :: Qref,fmin,fmax
double precision, dimension(1:N), intent(out) :: point,poids
integer k,m
logical, dimension(1:N) :: ia
double precision ref,freq,chi2
double precision, dimension(1:N,1:N) :: covar
double precision, dimension(1:2*N-1) :: x,y_ref,sig
m = 2*N-1
call remplit_point(fmin,fmax,N,point)
ref = 1.0 / Qref
do k=1,m
freq = (fmax/fmin) ** ((k - 1.)/(m - 1.))
freq = TWO_PI * fmin * freq
x(k) = freq
y_ref(k) = ref
sig(k) = 1.
enddo
do k=1,N
ia(k) = .TRUE.
enddo
call lfit_zener(x,y_ref,sig,m,poids,ia,covar,chi2,N,Qref,point)
END SUBROUTINE classical_linear_least_squares
! Calcul des coefficients par optimisation non-lineaire avec contraintes
SUBROUTINE solvopt(n,x,f,fun,flg,grad,options,flfc,func,flgc,gradc,Qref,Kopt,theta_min,theta_max,f_min,f_max)
!-----------------------------------------------------------------------------
! The subroutine SOLVOPT performs a modified version of Shor's r-algorithm in
! order to find a local minimum resp. maximum of a nonlinear function
! defined on the n-dimensional Euclidean space
! or
! a local minimum for a nonlinear constrained problem:
! min { f(x): g(x) (<)= 0, g(x) in R(m), x in R(n) }.
! Arguments:
! n is the space dimension (integer*4),
! x is the n-vector, the coordinates of the starting point
! at a call to the subroutine and the optimizer at regular return
! (double precision),
! f returns the optimum function value
! (double precision),
! fun is the entry name of a subroutine which computes the value
! of the function <fun> at a point x, should be declared as external
! in a calling routine,
! synopsis: fun(x,f)
! grad is the entry name of a subroutine which computes the gradient
! vector of the function <fun> at a point x, should be declared as
! external in a calling routine,
! synopsis: grad(x,g)
! func is the entry name of a subroutine which computes the MAXIMAL
! RESIDIAL!!! (a scalar) for a set of constraints at a point x,
! should be declared as external in a calling routine,
! synopsis: func(x,fc)
! gradc is the entry name of a subroutine which computes the gradient
! vector for a constraint with the MAXIMAL RESIDUAL at a point x,
! should be declared as external in a calling routine,
! synopsis: gradc(x,gc)
! flg, (logical) is a flag for the use of a subroutine <grad>:
! .true. means gradients are calculated by the user-supplied routine.
! flfc, (logical) is a flag for a constrained problem:
! .true. means the maximal residual for a set of constraints
! is calculated by <func>.
! flgc, (logical) is a flag for the use of a subroutine <gradc>:
! .true. means gradients of the constraints are calculated
! by the user-supplied routine.
! options is a vector of optional parameters (double precision):
! options(1)= H, where sign(H)=-1 resp. sign(H)=+1 means minimize resp.
! maximize <fun> (valid only for an unconstrained problem) and
! H itself is a factor for the initial trial step size
! (options(1)=-1.d0 by default),
! options(2)= relative error for the argument in terms of the infinity-norm
! (1.d-4 by default),
! options(3)= relative error for the function value (1.d-6 by default),
! options(4)= limit for the number of iterations (1.5d4 by default),
! options(5)= control of the display of intermediate results and error
! resp. warning messages (default value is 0.d0, i.e., no intermediate
! output but error and warning messages, see the manual for more),
! options(6)= maximal admissible residual for a set of constraints
! (options(6)=1.d-8 by default, see the manual for more),
! *options(7)= the coefficient of space dilation (2.5d0 by default),
! *options(8)= lower bound for the stepsize used for the difference
! approximation of gradients (1.d-11 by default,see the manual for more).
! (* ... changes should be done with care)
! Returned optional values:
! options(9), the number of iterations, if positive,
! or an abnormal stop code, if negative (see manual for more),
! -1: allocation error,
! -2: improper space dimension,
! -3: <fun> returns an improper value,
! -4: <grad> returns a zero vector or improper value at the
! starting point,
! -5: <func> returns an improper value,
! -6: <gradc> returns an improper value,
! -7: function is unbounded,
! -8: gradient is zero at the point,
! but stopping criteria are not fulfilled,
! -9: iterations limit exceeded,
! -11: Premature stop is possible,
! -12: Result may not provide the true optimum,
! -13: Function is flat: result may be inaccurate
! in view of a point.
! -14: Function is steep: result may be inaccurate
! in view of a function value,
! options(10), the number of objective function evaluations, and
! options(11), the number of gradient evaluations.
! options(12), the number of constraint function evaluations, and
! options(13), the number of constraint gradient evaluations.
! ____________________________________________________________________________
!
IMPLICIT NONE
!include 'messages.inc'
integer, intent(in) :: Kopt
double precision, intent(in) :: Qref,theta_min,theta_max,f_min,f_max
logical flg,flgc,flfc, constr, app, appconstr
logical FsbPnt, FsbPnt1, termflag, stopf
logical stopping, dispwarn, Reset, ksm,knan,obj
integer n, kstore, ajp,ajpp,knorms, k, kcheck, numelem
integer dispdata, ld, mxtc, termx, limxterm, nzero, krerun
integer warnno, kflat, stepvanish, i,j,ni,ii, kd,kj,kc,ip
integer iterlimit, kg,k1,k2, kless, allocerr
double precision options(13),doptions(13)
double precision x(n),f
double precision nsteps(3), gnorms(10), kk, nx
double precision ajb,ajs, des, dq,du20,du10,du03
double precision n_float, cnteps
double precision low_bound, ZeroGrad, ddx, y
double precision lowxbound, lowfbound, detfr, detxr, grbnd
double precision fp,fp1,fc,f1,f2,fm,fopt,frec,fst, fp_rate
double precision PenCoef, PenCoefNew
double precision gamma,w,wdef,h1,h,hp
double precision dx,ng,ngc,nng,ngt,nrmz,ng1,d,dd, laststep
double precision zero,one,two,three,four,five,six,seven
double precision eight,nine,ten,hundr
double precision infty, epsnorm,epsnorm2,powerm12
double precision, dimension(:,:), allocatable :: B
double precision, dimension(:), allocatable :: g
double precision, dimension(:), allocatable :: g0
double precision, dimension(:), allocatable :: g1
double precision, dimension(:), allocatable :: gt
double precision, dimension(:), allocatable :: gc
double precision, dimension(:), allocatable :: z
double precision, dimension(:), allocatable :: x1
double precision, dimension(:), allocatable :: xopt
double precision, dimension(:), allocatable :: xrec
double precision, dimension(:), allocatable :: grec
double precision, dimension(:), allocatable :: xx
double precision, dimension(:), allocatable :: deltax
integer, dimension(:), allocatable :: idx
character(len=100) :: endwarn
character(len=19) :: allocerrstr
external fun,grad,func,gradc
data zero/0.d0/, one/1.d0/, two/2.d0/, three/3.d0/, four/4.d0/, &
five/5.d0/, six/6.d0/, seven/7.d0/, eight/8.d0/, nine/9.d0/, &
ten/1.d1/, hundr/1.d2/, powerm12/1.d-12/, &
infty /1.d100/, epsnorm /1.d-15/, epsnorm2 /1.d-30/, &
allocerrstr/'Allocation Error = '/
! Check the dimension:
if (n<2) then
print *, 'SolvOpt error:'
print *, 'Improper space dimension.'
stop 'error in allocate statement in SolvOpt'
options(9)=-one
goto 999
endif
n_float=dble(n)
! allocate working arrays:
allocate (B(n,n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (g(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (g0(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (g1(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (gt(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (gc(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (z(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (x1(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (xopt(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (xrec(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (grec(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (xx(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (deltax(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
allocate (idx(n),stat=allocerr)
if (allocerr/=0) then
options(9)=-one
print *,allocerrstr,allocerr
stop 'error in allocate statement in SolvOpt'
endif
! store flags:
app=.not.flg
constr=flfc
appconstr=.not.flgc
! Default values for options:
call soptions(doptions)
do i=1,8
if (options(i)==zero) then
options(i)=doptions(i)
else if (i==2.or.i==3.or.i==6) then
options(i)=dmax1(options(i),powerm12)
options(i)=dmin1(options(i),one)
if (i==2)options(i)=dmax1(options(i),options(8)*hundr)
else if (i==7) then
options(7)=dmax1(options(i),1.5d0)
endif
enddo
! WORKING CONSTANTS AND COUNTERS ----{
options(10)=zero !! counter for function calculations
options(11)=zero !! counter for gradient calculations
options(12)=zero !! counter for constraint function calculations
options(13)=zero !! counter for constraint gradient calculations
iterlimit=idint(options(4))
if (constr) then
h1=-one !! NLP: restricted to minimization
cnteps=options(6)
else
h1=dsign(one,options(1)) !! Minimize resp. maximize a function
endif
k=0 !! Iteration counter
wdef=one/options(7)-one !! Default space transf. coeff.
! Gamma control ---{
ajb=one+1.d-1/n_float**2 !! Base I
ajp=20
ajpp=ajp !! Start value for the power
ajs=1.15d0 !! Base II
knorms=0
do i=1,10
gnorms(i)=zero
enddo
!---}
! Display control ---{
if (options(5)<=zero) then
dispdata=0
if (options(5)==-one) then
dispwarn=.false.
else
dispwarn=.true.
endif
else
dispdata=idnint(options(5))
dispwarn=.true.
endif
ld=dispdata
!---}
! Stepsize control ---{
dq=5.1d0 !! Step divider (at f_{i+1}>gamma*f_{i})
du20=two
du10=1.5d0
du03=1.05d0 !! Step multipliers (at certain steps made)
kstore=3
do i=1,kstore
nsteps(i)=zero !! Steps made at the last 'kstore' iterations
enddo
if (app) then
des=6.3d0 !! Desired number of steps per 1-D search
else
des=3.3d0
endif
mxtc=3 !! Number of trial cycles (steep wall detect)
!---}
termx=0
limxterm=50 !! Counter and limit for x-criterion
! stepsize for gradient approximation
ddx=dmax1(1.d-11,options(8))
low_bound=-one+1.d-4 !! Lower bound cosine used to detect a ravine
ZeroGrad=n_float*1.d-16 !! Lower bound for a gradient norm
nzero=0 !! Zero-gradient events counter
! Low bound for the values of variables to take into account
lowxbound=dmax1(options(2),1.d-3)
! Lower bound for function values to be considered as making difference
lowfbound=options(3)**2
krerun=0 !! Re-run events counter
detfr=options(3)*hundr !! Relative error for f/f_{record}
detxr=options(2)*ten !! Relative error for norm(x)/norm(x_{record})
warnno=0 !! the number of warn.mess. to end with
kflat=0 !! counter for points of flatness
stepvanish=0 !! counter for vanished steps
stopf=.false.
! ----} End of setting constants
! ----} End of the preamble
!--------------------------------------------------------------------
! COMPUTE THE FUNCTION ( FIRST TIME ) ----{
call fun(x,f,Qref,n/2,n,Kopt,f_min,f_max)
options(10)=options(10)+one
if (dabs(f)>=infty) then
if (dispwarn) then
print *,'SolvOpt error:'
print *,'Function equals infinity at the point.'
print *,'Choose another starting point.'
endif
options(9)=-three
goto 999
endif
do i=1,n
xrec(i)=x(i)
enddo
frec=f !! record point and function value
! Constrained problem
if (constr) then
kless=0
fp=f
call func(x,fc,n/2,n,theta_min,theta_max)
options(12)=options(12)+one
if (dabs(fc)>=infty) then
if (dispwarn) then
print *,'SolvOpt error:'
print *,'<FUNC> returns infinite value at the point.'
print *,'Choose another starting point.'
endif
options(9)=-five
goto 999
endif
PenCoef=one !! first rough approximation
if (fc<=cnteps) then
FsbPnt=.true. !! feasible point
fc=zero
else
FsbPnt=.false.
endif
f=f+PenCoef*fc
endif
! ----}
! COMPUTE THE GRADIENT ( FIRST TIME ) ----{
if (app) then
do i=1,n
deltax(i)=h1*ddx
enddo
obj=.true.
!if (constr) then
!call apprgrdn()
!else
!call apprgrdn()
!endif
options(10)=options(10)+n_float
else
call grad(x,g,Qref,n/2,n,Kopt,f_min,f_max)
options(11)=options(11)+one
endif
ng=zero
do i=1,n
ng=ng+g(i)*g(i)
enddo
ng=dsqrt(ng)
if (ng>=infty) then
if (dispwarn) then
print *,'SolvOpt error:'
print *,'Gradient equals infinity at the starting point.'
print *,'Choose another starting point.'
endif
options(9)=-four
goto 999
else if (ng<ZeroGrad) then
if (dispwarn) then
print *,'SolvOpt error:'
print *,'Gradient equals zero at the starting point.'
print *,'Choose another starting point.'
endif
options(9)=-four
goto 999
endif
if (constr) then
if (.not.FsbPnt) then
!if (appconstr) then
!do j=1,n
!if (x(j)>=zero) then
!deltax(j)=ddx
!else
!deltax(j)=-ddx
!endif
!enddo
!obj=.false.
!call apprgrdn()
if (.not. appconstr) then
call gradc(x,gc,n/2,n,theta_min,theta_max)
endif
ngc=zero
do i=1,n
ngc=ngc+gc(i)*gc(i)
enddo
ngc=dsqrt(ngc)
if (ng>=infty) then
if (dispwarn) then
print *,'SolvOpt error:'
print *,'<GRADC> returns infinite vector at the point.'
print *,'Choose another starting point.'
endif
options(9)=-six
goto 999
else if (ng<ZeroGrad) then
if (dispwarn) then
print *,'SolvOpt error:'
print *,'<GRADC> returns zero vector at an infeasible point.'
endif
options(9)=-six
goto 999
endif
do i=1,n
g(i)=g(i)+PenCoef*gc(i)
enddo
ng=zero
do i=1,n
ng=ng+g(i)*g(i)
grec(i)=g(i)
enddo
ng=dsqrt(ng)
endif
endif
do i=1,n
grec(i)=g(i)
enddo
nng=ng
! ----}
! INITIAL STEPSIZE
d=zero
do i=1,n
if (d<dabs(x(i))) d=dabs(x(i))
enddo
h=h1*dsqrt(options(2))*d !! smallest possible stepsize
if (dabs(options(1))/=one) then
h=h1*dmax1(dabs(options(1)),dabs(h)) !! user-supplied stepsize
else
h=h1*dmax1(one/dlog(ng+1.1d0),dabs(h)) !! calculated stepsize
endif
! RESETTING LOOP ----{
do while (.true.)
kcheck=0 !! Set checkpoint counter.
kg=0 !! stepsizes stored
kj=0 !! ravine jump counter
do i=1,n
do j=1,n
B(i,j)=zero
enddo
B(i,i)=one !! re-set transf. matrix to identity
g1(i)=g(i)
enddo
fst=f
dx=0
! ----}
! MAIN ITERATIONS ----{
do while (.true.)
k=k+1
kcheck=kcheck+1
laststep=dx
! ADJUST GAMMA --{
gamma=one+dmax1(ajb**((ajp-kcheck)*n),two*options(3))
gamma=dmin1 ( gamma,ajs**dmax1(one,dlog10(nng+one)) )
! --}
ngt=zero
ng1=zero
dd=zero
do i=1,n
d=zero
do j=1,n
d=d+B(j,i)*g(j)
enddo
gt(i)=d
dd=dd+d*g1(i)
ngt=ngt+d*d
ng1=ng1+g1(i)*g1(i)
enddo
ngt=dsqrt(ngt)
ng1=dsqrt(ng1)
dd=dd/ngt/ng1
w=wdef