-
Notifications
You must be signed in to change notification settings - Fork 0
/
minresqlpModule.f90
1382 lines (1261 loc) · 53.2 KB
/
minresqlpModule.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
!> \file
!! MINRES-QLP algorithm.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File minresqlpModule.f90
!
!> MINRESQLP solves symmetric systems Ax = b or min ||Ax - b||_2,
!! where the matrix A may be indefinite and/or singular.
!! "A" is really (A - shift*I), where shift is an input real scalar.
!!
!! \verbatim
!! 09 Sep 2013: Version 27
!!-------------------------------------------------------------------
!!
!! The software for MINRES-QLP is provided by SOL, Stanford University
!! under the terms of the OSI Common Public License (CPL)
!! http://www.opensource.org/licenses/cpl1.0.php
!! or the BSD License
!! http://www.opensource.org/licenses/bsd-license.php
!!
!!-------------------------------------------------------------------
!!
!! Authors:
!! Sou-Cheng Choi <[email protected]>
!! Computation Institute (CI)
!! University of Chicago
!! Chicago, IL 60637, USA
!!
!! Michael Saunders <[email protected]>
!! Systems Optimization Laboratory (SOL)
!! Stanford University
!! Stanford, CA 94305-4026, USA
!!
!! Contributor:
!! Christopher Paige <[email protected]>
!!
!! See also: Makefile
!! \endverbatim
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
module minresqlpModule
use minresqlpDataModule, only : dp, ip, one, zero, eps, realmin, prcsn, debug
use minresqlpBlasModule, only : dnrm2
implicit none
private ! sets default for module
public :: MINRESQLP, SYMORTHO
contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> Solution of linear equation system or least squares problem.
!!
!! \verbatim
!!------------------------------------------------------------------
!!
!! MINRESQLP is designed to solve the system of linear equations
!!
!! Ax = b
!!
!! or the least-squares problem
!!
!! min || Ax - b ||_2,
!!
!! where A is an n by n symmetric matrix and b is a given vector.
!! The matrix A may be indefinite and/or singular.
!!
!! 1. If A is known to be positive definite, the Conjugate Gradient
!! Method might be preferred, since it requires roughly the same
!! number of iterations as MINRESQLP but less work per iteration.
!! But if a low-accuracy solution is adequate, MINRESQLP will
!! terminate sooner.
!!
!! 2. If A is indefinite but Ax = b is known to have a solution
!! (e.g. if A is nonsingular), SYMMLQ might be preferred,
!! since it requires roughly the same number of iterations as
!! MINRESQLP but slightly less work per iteration.
!!
!! 3. If A is indefinite and well-conditioned, and Ax = b has a
!! solution, i.e., it is not a least-squares problem, MINRES might
!! be preferred as it requires the same number of iterations as
!! MINRESQLP but slightly less work per iteration.
!!
!! The matrix A is intended to be large and sparse. It is accessed
!! by means of a subroutine call of the form
!!
!! call Aprod ( n, x, y )
!!
!! which must return the product y = Ax for any given vector x.
!!
!!
!! More generally, MINRESQLP is designed to solve the system
!!
!! (A - shift*I) x = b
!! or
!! min || (A - shift*I) x - b ||_2,
!!
!! where shift is a specified real scalar. Again, the matrix
!! (A - shift*I) may be indefinite and/or singular.
!! The work per iteration is very slightly less if shift = 0.
!!
!! Note: If shift is an approximate eigenvalue of A
!! and b is an approximate eigenvector, x might prove to be
!! a better approximate eigenvector, as in the methods of
!! inverse iteration and/or Rayleigh-quotient iteration.
!! However, we are not yet sure on that -- it may be better
!! to use SYMMLQ.
!!
!! In this documentation, ' denotes the transpose of
!! a vector or a matrix.
!!
!! A further option is that of preconditioning, which may reduce
!! the number of iterations required. If M = C C' is a positive
!! definite matrix that is known to approximate (A - shift*I)
!! in some sense, and if systems of the form My = x can be
!! solved efficiently, the parameter Msolve may be used (see below).
!! When an external procedure Msolve is supplied, MINRESQLP will
!! implicitly solve the system of equations
!!
!! P (A - shift*I) P' xbar = P b,
!!
!! i.e. Abar xbar = bbar
!! where P = C**(-1),
!! Abar = P (A - shift*I) P',
!! bbar = P b,
!!
!! and return the solution x = P' xbar.
!! The associated residual is rbar = bbar - Abar xbar
!! = P (b - (A - shift*I)x)
!! = P r.
!!
!! In the discussion below, eps refers to the machine precision.
!! eps is computed by MINRESQLP. A typical value is eps = 2.22d-16
!! for IEEE double-precision arithmetic.
!!
!! Parameters
!! ----------
!! Some inputs are optional, with default values described below.
!! Mandatory inputs are n, Aprod, and b.
!! All outputs other than x are optional.
!!
!! n input The dimension of the matrix or operator A.
!!
!! b(n) input The rhs vector b.
!!
!! x(n) output Returns the computed solution x.
!!
!! Aprod external A subroutine defining the matrix A.
!! For a given vector x, the statement
!!
!! call Aprod ( n, x, y )
!!
!! must return the product y = Ax
!! without altering the vector x.
!! An extra call of Aprod is
!! used to check if A is symmetric.
!!
!! Msolve external An optional subroutine defining a
!! preconditioning matrix M, which should
!! approximate (A - shift*I) in some sense.
!! M must be positive definite.
!! For a given vector x, the statement
!!
!! call Msolve( n, x, y )
!!
!! must solve the linear system My = x
!! without altering the vector x.
!!
!! In general, M should be chosen so that Abar has
!! clustered eigenvalues. For example,
!! if A is positive definite, Abar would ideally
!! be close to a multiple of I.
!! If A or A - shift*I is indefinite, Abar might
!! be close to a multiple of diag( I -I ).
!!
!! shift input Should be zero if the system Ax = b is to be
!! solved. Otherwise, it could be an
!! approximation to an eigenvalue of A, such as
!! the Rayleigh quotient b'Ab / (b'b)
!! corresponding to the vector b.
!! If b is sufficiently like an eigenvector
!! corresponding to an eigenvalue near shift,
!! then the computed x may have very large
!! components. When normalized, x may be
!! closer to an eigenvector than b. Default to 0.
!!
!! nout input A file number. The calling program must open a file
!! for output using for example:
!! open(nout, file='MINRESQLP.txt', status='unknown')
!! If nout > 0, a summary of the iterations
!! will be printed on unit nout. If nout is absent or
!! the file associated with nout is not opened properly,
!! results will be written to 'MINRESQLP_tmp.txt'.
!! (Avoid 0, 5, 6 because by convention stderr=0,
!! stdin=5, stdout=6.)
!!
!! itnlim input An upper limit on the number of iterations. Default to 4n.
!!
!! rtol input A user-specified tolerance. MINRESQLP terminates
!! if it appears that norm(rbar) is smaller than
!! rtol*[norm(Abar)*norm(xbar) + norm(b)],
!! where rbar = bbar - Abar xbar,
!! or that norm(Abar*rbar) is smaller than
!! rtol*norm(Abar)*norm(rbar).
!!
!! If shift = 0 and Msolve is absent, MINRESQLP
!! terminates if norm(r) is smaller than
!! rtol*[norm(A)*norm(x) + norm(b)],
!! where r = b - Ax,
!! or if norm(A*r) is smaller than
!! rtol*norm(A)*norm(r).
!!
!! Default to machine precision.
!!
!! istop output An integer giving the reason for termination...
!! 0 Initial value of istop.
!!
!! 1 beta_{k+1} < eps.
!! Iteration k is the final Lanczos step.
!!
!! 2 beta2 = 0 in the Lanczos iteration; i.e. the
!! second Lanczos vector is zero. This means the
!! rhs is very special.
!!
!! If there is no preconditioner, b is an
!! eigenvector of Abar. Also, x = (1/alpha1) b
!! is a solution of Abar x = b.
!!
!! Otherwise (if Msolve is present), let My = b.
!! If shift is zero, y is a solution of the
!! generalized eigenvalue problem Ay = lambda My,
!! with lambda = alpha1 from the Lanczos vectors.
!!
!! In general, (A - shift*I)x = b
!! has the solution x = (1/alpha1) y
!! where My = b.
!!
!! 3 b = 0, so the exact solution is x = 0.
!! No iterations were performed.
!!
!! 4 Norm(rbar) appears to be less than
!! the value rtol * [ norm(Abar) * norm(xbar) + norm(b) ].
!! The solution in x should be an acceptable
!! solution of Abar x = b.
!!
!! 5 Norm(rbar) appears to be less than
!! the value eps * norm(Abar) * norm(xbar).
!! This means that the solution is as accurate as
!! seems reasonable on this machine.
!!
!! 6 Norm(Abar rbar) appears to be less than
!! the value rtol * norm(Abar) * norm(rbar).
!! The solution in x should be an acceptable
!! least-squares solution.
!!
!! 7 Norm(Abar rbar) appears to be less than
!! the value eps * norm(Abar) * norm(rbar).
!! This means that the least-squares solution is as
!! accurate as seems reasonable on this machine.
!!
!! 8 The iteration limit was reached before convergence.
!!
!! 9 The matrix defined by Aprod does not appear
!! to be symmetric.
!! For certain vectors y = Av and r = Ay, the
!! products y'y and r'v differ significantly.
!!
!! 10 The matrix defined by Msolve does not appear
!! to be symmetric.
!! For vectors satisfying My = v and Mr = y, the
!! products y'y and r'v differ significantly.
!!
!! 11 An inner product of the form x' M**(-1) x
!! was not positive, so the preconditioning matrix
!! M does not appear to be positive definite.
!!
!! 12 xnorm has exceeded maxxnorm or will exceed it
!! next iteration.
!!
!! 13 Acond (see below) has exceeded Acondlim or 0.1/eps,
!! so the matrix Abar must be very ill-conditioned.
!!
!! 14 | gamma_k | < eps.
!! This is very likely a least-squares problem but
!! x may not contain an acceptable solution yet.
!!
!! 15 norm(Abar x) < rtol * norm(Abar) * norm(x).
!! If disable = .true., then a null vector will be
!! obtained, given rtol.
!!
!! If istop >= 7, the final x may not be an
!! acceptable solution.
!!
!! itn output The number of iterations performed.
!!
!! Anorm output An estimate of the norm of the matrix operator
!! Abar = P (A - shift*I) P', where P = C**(-1).
!!
!! Acond output An estimate of the condition of Abar above.
!! This will usually be a substantial
!! under-estimate of the true condition.
!!
!! rnorm output An estimate of the norm of the final
!! transformed residual vector,
!! P (b - (A - shift*I) x).
!!
!! xnorm output An estimate of the norm of xbar.
!! This is sqrt( x'Mx ). If Msolve is absent,
!! xnorm is an estimate of norm(x).
!!
!! maxxnorm input An upper bound on norm(x). Default value is 1e7.
!!
!! trancond input If trancond > 1, a switch is made from MINRES
!! iterations to MINRES-QLP iterations when
!! Acond > trancond.
!! If trancond = 1, all iterations are MINRES-QLP
!! iterations.
!! If trancond = Acondlim, all iterations are
!! conventional MINRES iterations (which are
!! slightly cheaper).
!! Default to 1e7.
!!
!! Acondlim input An upper bound on Acond. Default value is 1e15.
!!
!! disable input All stopping conditions are disabled except
!! norm(Ax) / norm(x) < tol. Default to .false..
!!
!!------------------------------------------------------------------
!!
!! MINRESQLP is an implementation of the algorithm described in
!! the following references:
!!
!! Sou-Cheng Choi,
!! Iterative Methods for Singular Linear Equations and Least-
!! Squares Problems, PhD dissertation, ICME, Stanford University,
!! 2006.
!!
!! Sou-Cheng Choi, Christopher Paige, and Michael Saunders,
!! MINRES-QLP: A Krylov subspace method for indefinite or
!! singular symmetric systems, SIAM Journal of Scientific
!! Computing 33:4 (2011) 1810-1836.
!!
!! Sou-Cheng Choi and Michael Saunders,
!! ALGORITHM & DOCUMENTATION: MINRES-QLP for singular Symmetric and Hermitian
!! linear equations and least-squares problems, Technical Report,
!! ANL/MCS-P3027-0812, Computation Institute,
!! University of Chicago/Argonne National Laboratory, 2012.
!!
!! Sou-Cheng Choi and Michael Saunders,
!! ALGORITHM xxx: MINRES-QLP for singular Symmetric and Hermitian
!! linear equations and least-squares problems,
!! ACM Transactions on Mathematical Software, to appear, 2013.
!!
!! FORTRAN 90 and MATLAB implementations are
!! downloadable from
!! http://www.stanford.edu/group/SOL/software.html
!! http://home.uchicago.edu/sctchoi/
!!
!!------------------------------------------------------------------
!!
!! MINRESQLP development:
!! 14 Dec 2006: Sou-Cheng's thesis completed.
!! MINRESQLP includes a stopping rule for singular
!! systems (using an estimate of ||Ar||) and very many
!! other things(!).
!! Note that ||Ar|| small => r is a null vector for A.
!! 09 Oct 2007: F90 version constructed from the F77 version.
!! Initially used compiler option -r8, but this is
!! nonstandard.
!! 15 Oct 2007: Test on Arnorm = ||Ar|| added to recognize
!! singular systems.
!! 15 Oct 2007: Temporarily used real(8) everywhere.
!! 16 Oct 2007: Use minresqlpDataModule to define
!! dp = selected_real_kind(15).
!! We need "use minresqlpDataModule" at the
!! beginning of modules AND inside interfaces.
!! 06 Jun 2010: Added comments.
!! 12 Jul 2011: Created complex version zminresqlpModule.f90
!! from real version minresqlpModule.f90.
!! 23 Aug 2011: (1) Tim Hopkins ran version 17 on the NAG Fortran compiler
!! We removed half a dozen unused variables in MINRESQLP
!! and also local var sgn_a and sgn_b in SMMORTHO,
!! as they result in division by zero for inputs a=b=0.
!! (2) Version 18 was submitted to ACM TOMS for review.
!! 20 Aug 2012: Version 19:
!! (1) Added optional inputs and outputs, and
!! default values for optional inputs.
!! (2) Removed inputs 'checkA' and 'precon'.
!! (3) Changed slightly the order of parameters in the
!! MINRESQLP API.
!! (4) Updated documentation.
!! (5) Fixed a minor bug in printing x(1) in iteration
!! log during MINRES mode.
!! (6) Made sure MINRESQLP is portable in both single
!! and double precison.
!! (7) Fixed a bug to ensure the 2x2 Hermitian reflectors
!! are orthonormal. Make output c real.
!! 24 Apr 2013: istop = 12 now means xnorm just exceeded maxxnorm.
!! 28 Jun 2013: likeLS introduced to terminate with big xnorm
!! only if the problem seems to be singular and inconsistent.
!! 08 Jul 2013: (1) dot_product replaces ddotc.
!! 04 Aug 2013: If present(maxxnorm), use maxxnorm_ = min(maxxnorm, one/eps).
!! 09 Sep 2013: Initialize relresl and relAresl to zero.
!!------------------------------------------------------------------
!! \endverbatim
subroutine MINRESQLP( n, Aprod, b, shift, Msolve, disable, nout, &
itnlim, rtol, maxxnorm, trancond, Acondlim, &
x, istop, itn, rnorm, Arnorm, xnorm, Anorm, Acond )
! Inputs
integer(ip), intent(in) :: n
real(dp), intent(in) :: b(n)
integer(ip), intent(in), optional :: itnlim, nout
logical, intent(in), optional :: disable
real(dp), intent(in), optional :: shift
real(dp), intent(in), optional :: rtol, maxxnorm, trancond, Acondlim
! Outputs
real(dp), intent(out) :: x(n)
integer(ip), intent(out), optional :: istop, itn
real(dp), intent(out), optional :: rnorm, Arnorm, xnorm, Anorm, Acond
optional :: Msolve
interface
subroutine Aprod(n,x,y) ! y := A*x
use minresqlpDataModule
integer(ip), intent(in) :: n
real(dp), intent(in) :: x(n)
real(dp), intent(out) :: y(n)
end subroutine Aprod
subroutine Msolve(n,x,y) ! Solve M*y = x
use minresqlpDataModule
integer(ip), intent(in) :: n
real(dp), intent(in) :: x(n)
real(dp), intent(out) :: y(n)
end subroutine Msolve
end interface
intrinsic :: abs, sqrt, present, floor, log10, dot_product
! Local arrays and variables
real(dp) :: shift_
real(dp) :: rtol_, maxxnorm_, trancond_, Acondlim_
real(dp) :: rnorm_, Arnorm_, xnorm_, Anorm_, Acond_
logical :: checkA_, precon_, disable_
integer(ip) :: itnlim_, nout_, istop_, itn_
real(dp) :: r1(n), r2(n), v(n), w(n), wl(n), wl2(n),&
xl2(n), y(n), vec2(2), vec3(3)
real(dp) :: abs_gama, Acondl, alfa, Anorml, Arnorml,&
Axnorm, beta, beta1, betal, betan, &
epsa, epsx, gminl2, ieps, pnorm, &
relAres, relAresl, relres, relresl, &
rnorml, rootl, t1, t2, xl2norm, &
xnorm_tmp, xnorml, z, &
cr1, cr2, cs, dbar, dlta, dlta_QLP, &
dlta_tmp, dltan, epln, eplnn, eta, etal,&
etal2, gama, gama_QLP, gama_tmp, gamal, &
gamal_QLP, gamal_tmp, gamal2, gamal3, &
gbar, gmin, gminl, phi, s, sn, sr1, sr2,&
t, tau, taul, taul2, u, u_QLP, &
ul, ul_QLP, ul2, ul2_QLP, ul3, ul4, &
vepln, vepln_QLP, veplnl, veplnl2, x1last
integer(ip) :: j, QLPiter, headlines, lines, nprint, flag0, ios
logical :: prnt, done, lastiter, connected, named_file, likeLS
character(len=20) :: filename
character(len=2) :: QLPstr = ' '
! Local constants
real(dp), parameter :: EPSINV = 10.0_dp**floor(log10(one/eps))
real(dp) :: NORMMAX = 10.0_dp**floor(log10(one/eps)/2)
character(len=*), parameter :: enter = ' Enter MINRES-QLP. '
character(len=*), parameter :: exitt = ' Exit MINRES-QLP. '
character(len=*), parameter :: msg(1:15) = &
(/ 'beta_{k+1} < eps. ', & ! 1
'beta2 = 0. If M = I, b and x are eigenvectors of A. ', & ! 2
'beta1 = 0. The exact solution is x = 0. ', & ! 3
'A solution to (poss. singular) Ax = b found, given rtol. ', & ! 4
'A solution to (poss. singular) Ax = b found, given eps. ', & ! 5
'Pseudoinverse solution for singular LS problem, given rtol. ', & ! 6
'Pseudoinverse solution for singular LS problem, given eps. ', & ! 7
'The iteration limit was reached. ', & ! 8
'The operator defined by Aprod appears to be unsymmetric. ', & ! 9
'The operator defined by Msolve appears to be unsymmetric. ', & ! 10
'The operator defined by Msolve appears to be indefinite. ', & ! 11
'xnorm has exceeded maxxnorm or will exceed it next iteration. ', & ! 12
'Acond has exceeded Acondlim or 0.1/eps. ', & ! 13
'Least-squares problem but no converged solution yet. ', & ! 14
'A null vector obtained, given rtol. ' /) ! 15
character(len=*), parameter :: ddebugStr1 = "(a, T5, i0, a, 5(e12.3))"
character(len=*), parameter :: ddebugStr2 = "(5(a, i0, a, e12.3, a))"
character(len=*), parameter :: headerStr = &
"(// 1p, a, 4x, 'Solution of symmetric Ax = b'" // &
" / ' n =', i7, 6x, '||b|| =', e11.2, 3x," // &
" 'precon =', l4 " // &
" / ' itnlim =', i7, 6x, 'rtol =', e11.2, 3x," // &
" 'shift =', e23.14 " // &
" / ' maxxnorm =', e11.2, 2x, 'Acondlim =', e11.2, 3x,"// &
" 'trancond =', e11.2)"
character(len=*), parameter :: tableHeaderStr = &
"(// ' iter x(1) xnorm rnorm Arnorm '," // &
" 'Compatible LS norm(A) cond(A)')"
character(len=*), parameter :: itnStr = "(1p, i8, e19.10, 7e10.2, a)"
character(len=*), parameter :: finalStr1 = &
"(/ 1p, a, 5x, a, i3, 14x, a, i8 " // &
" / a, 5x, a, e12.4, 5x, a, e12.4 " // &
" / a, 5x, a, e12.4, 5x, a, e12.4 " // &
" / a, 5x, a, e12.4 )"
character(len=*), parameter :: finalStr2 = "( a, 5x, a )"
!------------------------------------------------------------------
prnt = .FALSE.
t1 = zero
t2 = zero
gminl = zero
! Optional inputs
if (present(shift)) then
shift_ = shift
else
shift_ = zero
end if
checkA_ = .true.
if (present(disable)) then
disable_ = disable
else
disable_ = .false.
end if
if (present(itnlim)) then
itnlim_ = itnlim
else
itnlim_ = 4 * n
end if
connected = .false.
filename = "MINRESQLP_tmp.txt"
nout_ = 10
if (present(nout)) then
nout_ = nout
inquire(unit=nout, opened=connected, named=named_file, name=filename)
!write(*,*) connected, named_file, filename
if (.not. connected) then
write(*,*) "File unit 'nout' is not open."
if (nout==5 .or. nout == 6) then
nout_ = 10
end if
end if
end if
if (.not. connected) then
write(*,*) 'nout_ = ', nout_
open(nout_, file=filename, status='unknown', iostat=ios)
write(*,*) 'ios = ', ios
if (ios /= 0) then
write(*,*) "Error opening file '", filename, "'."
STOP
end if
end if
if (present(rtol)) then
rtol_ = rtol
else
rtol_ = eps
end if
if (prcsn == 6) then
NORMMAX = 1.0e4_dp
end if
if (present(maxxnorm)) then
maxxnorm_ = min(maxxnorm, one/eps)
else
maxxnorm_ = NORMMAX
end if
if (present(trancond)) then
trancond_ = min(trancond, NORMMAX)
else
trancond_ = NORMMAX
end if
if (present(Acondlim)) then
Acondlim_ = min(Acondlim, EPSINV)
else
Acondlim_ = EPSINV
end if
if (present(Msolve)) then
precon_ = .true.
else
precon_ = .false.
end if
!------------------------------------------------------------------
! Print heading and initialize.
!------------------------------------------------------------------
nprint = min(n,20)
!debug = .true.
lastiter = .false.
flag0 = 0
istop_ = flag0
beta1 = dnrm2(n, b, 1)
ieps = 0.1_dp/eps
itn_ = 0
QLPiter = 0
xnorm_ = zero
xl2norm = zero
Axnorm = zero
Anorm_ = zero
Acond_ = one
pnorm = zero
relresl = zero
relAresl = zero
x = zero
xl2 = zero
x1last = x(1)
if (nout_ > 0) then
write(nout_, headerStr) enter, n, beta1, precon_, itnlim_, rtol_, &
shift_, maxxnorm_, Acondlim_, trancond_
end if
!------------------------------------------------------------------
! Set up y and v for the first Lanczos vector v1.
! y = beta1 P'v1, where P = C**(-1).
! v is really P'v1.
!------------------------------------------------------------------
y = b
r1 = b
if ( precon_ ) then
call Msolve( n, b, y )
end if
beta1 = dot_product(b, y)
if (beta1 < zero .and. dnrm2(n, y, 1) > eps) then ! M must be indefinite.
istop_ = 11
end if
if (beta1 == zero) then ! b = 0 exactly. Stop with x = 0.
istop_ = 3
end if
beta1 = sqrt( beta1 ) ! Normalize y to get v1 later.
if (debug) then
write(*,ddebugStr1) ' y_', itn_, ' = ', (y(j), j=1,nprint)
write(*,*) 'beta1 ', beta1
end if
!------------------------------------------------------------------
! See if Msolve is symmetric.
!------------------------------------------------------------------
if (checkA_ .and. precon_) then
call Msolve( n, y, r2 )
s = dot_product( y, y )
t = dot_product(r1, r2)
z = abs( s - t )
epsa = (abs(s) + eps) * eps**0.33333_dp
if (z > epsa) then
istop_ = 10
end if
end if
!------------------------------------------------------------------
! See if Aprod is symmetric.
!------------------------------------------------------------------
if (checkA_) then
call Aprod ( n, y, w ) ! w = A*y
call Aprod ( n, w, r2 ) ! r2 = A*w
s = dot_product( w, w )
t = dot_product( y, r2)
z = abs( s - t )
epsa = (abs(s) + eps) * eps**0.33333_dp
if (z > epsa) then
istop_ = 9
end if
end if
!------------------------------------------------------------------
! Initialize other quantities.
!------------------------------------------------------------------
tau = zero
taul = zero
gmin = zero
beta = zero
betan = beta1
dbar = zero
dltan = zero
eta = zero
etal = zero
etal2 = zero
vepln = zero
veplnl = zero
veplnl2= zero
eplnn = zero
gama = zero
gamal = zero
gamal2 = zero
phi = beta1
cr1 = -one
sr1 = zero
cr2 = -one
sr2 = zero
cs = -one
sn = zero
ul3 = zero
ul2 = zero
ul = zero
u = zero
lines = 1
headlines = 30 * lines
rnorm_ = betan
relres = rnorm_ / (Anorm_*xnorm_ + beta1)
relAres= zero
r2 = b
w = zero ! vector of zeros
wl = zero
done = .false.
if (debug) then
write(*,*)
write(*,*) 'Checking variable values before main loop'
write(*,*) 'istop ', istop_, ' done ', done, ' itn ', itn_, ' QLPiter' , QLPiter
write(*,*) 'lastiter', lastiter, ' lines ', lines, ' headlines ', headlines
write(*,*) 'beta ', beta, ' tau ', tau, ' taul ', taul, ' phi ', phi
write(*,*) 'betan ', betan, ' gmin ', gmin, ' cs ', cs, ' sn ', sn
write(*,*) 'cr1 ', cr1, ' sr1 ', sr1, ' cr2 ', cr2, ' sr2 ', sr2
write(*,*) 'dltan ', dltan, ' eplnn ', eplnn, ' gama ', gama, ' gamal ', gamal
write(*,*) 'gamal2 ', gamal2, ' eta ', eta, ' etal ', etal, 'etal2', etal2
write(*,*) 'vepln ', vepln, ' veplnl', veplnl, ' veplnl2 ', veplnl2, ' ul3 ', ul3
write(*,*) 'ul2 ', ul2, ' ul ', ul, ' u ', u, ' rnorm ', rnorm_
write(*,*) 'xnorm ', xnorm_, ' xl2norm ', xl2norm, ' pnorm ', pnorm, ' Axnorm ', Axnorm
write(*,*) 'Anorm ', Anorm_, ' Acond ', Acond_, ' relres ', relres
write(*,ddebugStr1) 'w_', itn_-1, ' = ', (wl(j), j=1,nprint)
write(*,ddebugStr1) 'w_', itn_, ' = ', (w(j), j=1,nprint)
write(*,ddebugStr1) 'x_', itn_, ' = ', (x(j), j=1,nprint)
end if
!------------------------------------------------------------------
! Main iteration loop.
!------------------------------------------------------------------
do while (istop_ <= flag0)
itn_ = itn_ + 1 ! k = itn = 1 first time through
!-----------------------------------------------------------------
! Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
! The general iteration is similar to the case k = 1 with v0 = 0:
!
! p1 = Operator * v1 - beta1 * v0,
! alpha1 = v1'p1,
! q2 = p2 - alpha1 * v1,
! beta2^2 = q2'q2,
! v2 = (1/beta2) q2.
!
! Again, y = betak P vk, where P = C**(-1).
! .... more description needed.
!-----------------------------------------------------------------
betal = beta; ! betal = betak
beta = betan;
s = one / beta ! Normalize previous vector (in y).
v = s*y; ! v = vk if P = I.
call Aprod ( n, v, y )
if (shift_ /= zero) then
y = y - shift_ * v
end if
if (itn_ >= 2) then
y = y + (- beta/betal) * r1
end if
alfa = dot_product(v, y) ! alphak
y = y + (- alfa/beta) * r2
r1 = r2
r2 = y
if ( .not. precon_ ) then
betan = dnrm2(n, y, 1) ! betan = ||y||_2
else
call Msolve( n, r2, y )
betan = dot_product(r2, y) ! betan = betak+1^2
if (betan > zero) then
betan = sqrt(betan)
elseif ( dnrm2(n, y, 1) > eps ) then ! M must be indefinite.
istop_ = 11
exit
end if
end if
if (itn_ == 1) then
vec2(1) = alfa
vec2(2) = betan
pnorm = dnrm2(2, vec2, 1)
else
vec3(1) = beta
vec3(2) = alfa
vec3(3) = betan
pnorm = dnrm2(3, vec3, 1)
end if
if (debug) then
write(*,*)
write(*,*) 'Lanczos iteration ', itn_
write(*,ddebugStr1) 'v_', itn_, ' = ', (v(j), j=1,nprint)
write(*,ddebugStr1) 'r1_', itn_, ' = ', (r1(j), j=1,nprint)
write(*,ddebugStr1) 'r2_', itn_, ' = ', (r2(j), j=1,nprint)
write(*,ddebugStr1) 'y_', itn_, ' = ', (y(j), j=1,nprint)
write(*,ddebugStr2) 'alpha_', itn_, ' = ', alfa, ', ', &
'beta_', itn_, ' = ', beta, ', ', &
'beta_', itn_+1, ' = ', betan, ', ', &
'pnorm_', itn_, ' = ', pnorm
end if
! Apply previous left reflection Qk-1 to get
! [deltak epslnk+1] = [cs sn][dbark 0 ]
! [gbar k dbar k+1] [sn -cs][alfak betak+1].
dbar = dltan
dlta = cs * dbar + sn * alfa ! dlta1 = 0 deltak
epln = eplnn;
gbar = sn * dbar - cs * alfa ! gbar 1 = alfa1 gbar k
eplnn = sn * betan ! eplnn2 = 0 epslnk+1
dltan = - cs * betan ! dbar 2 = beta2 dbar k+1
dlta_QLP = dlta;
if (debug) then
write(*,*)
write(*,*) 'Apply previous left reflection Q_{', itn_-1, ',', itn_, '}'
write(*,ddebugStr2) 'c_', itn_-1, ' = ', cs, ', ', &
's_', itn_-1, ' = ', sn
write(*,ddebugStr2) 'dlta_', itn_, ' = ', dlta, ', ', &
'gbar_', itn_, ' = ', gbar, ', ', &
'epln_', itn_+1, ' = ', eplnn, ', ', &
'dbar_', itn_+1, ' = ', dltan
end if
! Compute the current left reflection Qk
gamal3 = gamal2
gamal2 = gamal
gamal = gama
call symortho(gbar, betan, cs, sn, gama)
gama_tmp = gama;
taul2 = taul
taul = tau
tau = cs * phi
phi = sn * phi ! phik
Axnorm = sqrt( Axnorm**2 + tau**2 );
if (debug) then
write(*,*)
write(*,*) 'Compute the current left reflection Q_{', itn_, ',', itn_+1, '}'
write(*,ddebugStr2) 'c_', itn_, ' = ', cs, ', ', &
's_', itn_, ' = ', sn
write(*,ddebugStr2) 'tau_', itn_, ' = ', tau, ', ', &
'phi_', itn_, ' = ', phi, ', ', &
'gama_', itn_, ' = ', gama
end if
! Apply the previous right reflection P{k-2,k}
if (itn_ > 2) then
veplnl2 = veplnl
etal2 = etal
etal = eta
dlta_tmp = sr2 * vepln - cr2 * dlta
veplnl = cr2 * vepln + sr2 * dlta
dlta = dlta_tmp
eta = sr2 * gama
gama = -cr2 * gama
if (debug) then
write(*,*)
write(*,*) 'Apply the previous right reflection P_{', itn_-2, ',', itn_, '}'
write(*,ddebugStr2) 'cr2_', itn_, ' = ', cr2, ', ', &
'sr2_', itn_, ' = ', sr2
write(*,ddebugStr2) 'gamal2_', itn_, ' = ', gamal2, ', ', &
'gamal_', itn_, ' = ', gamal, ', ', &
'gama_', itn_, ' = ', gama
write(*,ddebugStr2) 'dlta_', itn_, ' = ', dlta, ', ', &
'vepln_', itn_-1, ' = ', veplnl
end if
end if
! Compute the current right reflection P{k-1,k}, P_12, P_23,...
if (itn_ > 1) then
call symortho(gamal, dlta, cr1, sr1, gamal_tmp)
gamal = gamal_tmp
vepln = sr1 * gama
gama = -cr1 * gama
if (debug) then
write(*,*)
write(*,*) 'Apply the current right reflection P_{', itn_-1, ',', itn_, '}'
write(*,ddebugStr2) 'cr1_ ', itn_, ' = ', cr1, ', ', &
'sr1_' , itn_, ' = ', sr1
write(*,ddebugStr2) 'gama_', itn_-2, ' = ', gamal2, ', ', &
'gama_', itn_-1, ' = ', gamal, ', ', &
'gama_', itn_, ' = ', gama
write(*,ddebugStr2) 'dlta_', itn_, ' = ', dlta, ', ', &
'vepln_', itn_-1, ' = ', veplnl, ', ', &
'eta_', itn_, ' = ', eta
end if
end if
! Update xnorm
xnorml = xnorm_
ul4 = ul3
ul3 = ul2
if (itn_ > 2) then
ul2 = ( taul2 - etal2 * ul4 - veplnl2 * ul3 ) / gamal2
if (debug) then
write(*,ddebugStr2) 'tau_', itn_-2, ' = ', taul2, ', ', &
'eta_', itn_-2, ' = ', etal2, ', ', &
'vepln_', itn_-2, ' = ', veplnl2, ', ', &
'gama_', itn_-2, ' = ', gamal2
end if
end if
if (itn_ > 1) then
ul = ( taul - etal * ul3 - veplnl * ul2) / gamal
if (debug) then
write(*,ddebugStr2) 'tau_', itn_-1, ' = ', taul, ', ', &
'eta_', itn_-1, ' = ', etal, ', ', &
'vepln_', itn_-1, ' = ', veplnl, ', ', &
'gamal_', itn_-1, ' = ', gamal
end if
end if
vec3(1) = xl2norm
vec3(2) = ul2
vec3(3) = ul
xnorm_tmp = dnrm2(3, vec3, 1) ! norm([xl2norm ul2 ul]);
if (abs(gama) > eps) then ! .and. xnorm_tmp < maxxnorm_) then
if (debug) then
write(*,ddebugStr2) 'tau_', itn_, ' = ', tau, ', ', &
'eta_', itn_, ' = ', eta, ', ', &
'vepln_', itn_, ' = ', vepln, ', ', &
'gama_', itn_, ' = ', gama
end if
u = (tau - eta*ul2 - vepln*ul) / gama
likeLS = relAresl < relresl
vec2(1) = xnorm_tmp
vec2(2) = u
if (likeLS .and. dnrm2(2, vec2, 1) > maxxnorm_) then
u = zero
istop_ = 12
end if
else
u = zero
istop_ = 14
end if
vec2(1) = xl2norm
vec2(2) = ul2
xl2norm = dnrm2(2, vec2, 1)
vec3(1) = xl2norm
vec3(2) = ul
vec3(3) = u
xnorm_ = dnrm2(3, vec3, 1)
if (Acond_ < trancond_ .and. istop_ == flag0 .and. QLPiter == 0) then !! MINRES updates
wl2 = wl
wl = w
if (gama_tmp > eps) then
s = one / gama_tmp
w = (v - epln*wl2 - dlta_QLP*wl) * s
end if
if (xnorm_ < maxxnorm_) then
x1last = x(1)
x = x + tau*w
else
istop_ = 12
lastiter = .true.
end if
if (debug) then
write(*,*)
write(*,*) 'MINRES updates'
write(*,ddebugStr2) 'gama_tmp_', itn_, ' = ', gama_tmp, ', ', &
'tau_', itn_, ' = ', tau, ', ', &
'epln_', itn_, ' = ', epln, ', ', &
'dlta_QLP_', itn_, ' = ', dlta_QLP
write(*,ddebugStr1) 'v_', itn_ , ' = ', (v(j), j=1,nprint)
write(*,ddebugStr1) 'w_', itn_ , ' = ', (w(j), j=1,nprint)
end if
else !! MINRES-QLP updates
QLPiter = QLPiter + 1;