Millepede-II V04-17-03
minresqlpModule.f90
Go to the documentation of this file.
1
3
4
5!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6! File minresqlpModule.f90
7!
40!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
41
43
45 use minresqlpblasmodule, only : dnrm2
46
47 implicit none
48
49 private ! sets default for module
50 public :: minresqlp, symortho
51
52contains
53
54 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
55
410
411 subroutine minresqlp( n, Aprod, b, shift, Msolve, disable, nout, &
412 itnlim, rtol, maxxnorm, trancond, Acondlim, &
413 x, istop, itn, rnorm, Arnorm, xnorm, Anorm, Acond )
414 ! Inputs
415 integer(ip), intent(in) :: n
416 real(dp), intent(in) :: b(n)
417 integer(ip), intent(in), optional :: itnlim, nout
418 logical, intent(in), optional :: disable
419 real(dp), intent(in), optional :: shift
420 real(dp), intent(in), optional :: rtol, maxxnorm, trancond, acondlim
421
422 ! Outputs
423 real(dp), intent(out) :: x(n)
424 integer(ip), intent(out), optional :: istop, itn
425 real(dp), intent(out), optional :: rnorm, arnorm, xnorm, anorm, acond
426
427 optional :: msolve
428
429 interface
430 subroutine aprod(n,x,y) ! y := A*x
432 integer(ip), intent(in) :: n
433 real(dp), intent(in) :: x(n)
434 real(dp), intent(out) :: y(n)
435 end subroutine aprod
436
437 subroutine msolve(n,x,y) ! Solve M*y = x
439 integer(ip), intent(in) :: n
440 real(dp), intent(in) :: x(n)
441 real(dp), intent(out) :: y(n)
442 end subroutine msolve
443 end interface
444
445 !obsolete intrinsic :: abs, sqrt, present, floor, log10, dot_product
446
447 ! Local arrays and variables
448 real(dp) :: shift_
449 real(dp) :: rtol_, maxxnorm_, trancond_, acondlim_
450 real(dp) :: rnorm_, arnorm_, xnorm_, anorm_, acond_
451 logical :: checka_, precon_, disable_
452 integer(ip) :: itnlim_, nout_, istop_, itn_
453
454 real(dp) :: r1(n), r2(n), v(n), w(n), wl(n), wl2(n),&
455 xl2(n), y(n), vec2(2), vec3(3)
456 real(dp) :: abs_gama, acondl, alfa, anorml, arnorml,&
457 axnorm, beta, beta1, betal, betan, &
458 epsa, epsx, gminl2, ieps, pnorm, &
459 relares, relaresl, relres, relresl, &
460 rnorml, rootl, t1, t2, xl2norm, &
461 xnorm_tmp, xnorml, z, &
462 cr1, cr2, cs, dbar, dlta, dlta_qlp, &
463 dlta_tmp, dltan, epln, eplnn, eta, etal,&
464 etal2, gama, gama_qlp, gama_tmp, gamal, &
465 gamal_qlp, gamal_tmp, gamal2, gamal3, &
466 gbar, gmin, gminl, phi, s, sn, sr1, sr2,&
467 t, tau, taul, taul2, u, u_qlp, &
468 ul, ul_qlp, ul2, ul2_qlp, ul3, ul4, &
469 vepln, vepln_qlp, veplnl, veplnl2, x1last
470
471 integer(ip) :: j, qlpiter, headlines, lines, nprint, flag0, ios
472 logical :: prnt, done, lastiter, connected, named_file, likels
473 character(len=20) :: filename
474 character(len=2) :: qlpstr = ' '
475
476 ! Local constants
477 real(dp), parameter :: epsinv = 10.0_dp**floor(log10(one/eps))
478 real(dp) :: normmax = 10.0_dp**floor(log10(one/eps)/2)
479 character(len=*), parameter :: enter = ' Enter MINRES-QLP. '
480 character(len=*), parameter :: exitt = ' Exit MINRES-QLP. '
481 character(len=*), parameter :: msg(1:15) = &
482 (/ 'beta_{k+1} < eps. ', & ! 1
483 'beta2 = 0. If M = I, b and x are eigenvectors of A. ', & ! 2
484 'beta1 = 0. The exact solution is x = 0. ', & ! 3
485 'A solution to (poss. singular) Ax = b found, given rtol. ', & ! 4
486 'A solution to (poss. singular) Ax = b found, given eps. ', & ! 5
487 'Pseudoinverse solution for singular LS problem, given rtol. ', & ! 6
488 'Pseudoinverse solution for singular LS problem, given eps. ', & ! 7
489 'The iteration limit was reached. ', & ! 8
490 'The operator defined by Aprod appears to be unsymmetric. ', & ! 9
491 'The operator defined by Msolve appears to be unsymmetric. ', & ! 10
492 'The operator defined by Msolve appears to be indefinite. ', & ! 11
493 'xnorm has exceeded maxxnorm or will exceed it next iteration. ', & ! 12
494 'Acond has exceeded Acondlim or 0.1/eps. ', & ! 13
495 'Least-squares problem but no converged solution yet. ', & ! 14
496 'A null vector obtained, given rtol. ' /) ! 15
497
498 character(len=*), parameter :: ddebugstr1 = "(a, T5, i0, a, 5(e12.3))"
499 character(len=*), parameter :: ddebugstr2 = "(5(a, i0, a, e12.3, a))"
500
501 character(len=*), parameter :: headerstr = &
502 "(// 1p, a, 4x, 'Solution of symmetric Ax = b'" // &
503 " / ' n =', i7, 6x, '||b|| =', e11.2, 3x," // &
504 " 'precon =', l4 " // &
505 " / ' itnlim =', i7, 6x, 'rtol =', e11.2, 3x," // &
506 " 'shift =', e23.14 " // &
507 " / ' maxxnorm =', e11.2, 2x, 'Acondlim =', e11.2, 3x,"// &
508 " 'trancond =', e11.2)"
509 character(len=*), parameter :: tableheaderstr = &
510 "(// ' iter x(1) xnorm rnorm Arnorm '," // &
511 " 'Compatible LS norm(A) cond(A)')"
512 character(len=*), parameter :: itnstr = "(1p, i8, e19.10, 7e10.2, a)"
513 character(len=*), parameter :: finalstr1 = &
514 "(/ 1p, a, 5x, a, i3, 14x, a, i8 " // &
515 " / a, 5x, a, e12.4, 5x, a, e12.4 " // &
516 " / a, 5x, a, e12.4, 5x, a, e12.4 " // &
517 " / a, 5x, a, e12.4 )"
518 character(len=*), parameter :: finalstr2 = "( a, 5x, a )"
519
520 !------------------------------------------------------------------
521 prnt = .false.
522 t1 = zero
523 t2 = zero
524 gminl = zero
525 ! Optional inputs
526 if (present(shift)) then
527 shift_ = shift
528 else
529 shift_ = zero
530 end if
531
532 checka_ = .true.
533
534 if (present(disable)) then
535 disable_ = disable
536 else
537 disable_ = .false.
538 end if
539
540 if (present(itnlim)) then
541 itnlim_ = itnlim
542 else
543 itnlim_ = 4 * n
544 end if
545
546 connected = .false.
547 filename = "MINRESQLP_tmp.txt"
548 nout_ = 10
549
550 if (present(nout)) then
551 nout_ = nout
552 inquire(unit=nout, opened=connected, named=named_file, name=filename)
553 !write(*,*) connected, named_file, filename
554 if (.not. connected) then
555 write(*,*) "File unit 'nout' is not open."
556 if (nout==5 .or. nout == 6) then
557 nout_ = 10
558 end if
559 end if
560 end if
561
562 if (.not. connected) then
563 write(*,*) 'nout_ = ', nout_
564 open(nout_, file=filename, status='unknown', iostat=ios)
565 write(*,*) 'ios = ', ios
566 if (ios /= 0) then
567 write(*,*) "Error opening file '", filename, "'."
568 stop
569 end if
570 end if
571
572 if (present(rtol)) then
573 rtol_ = rtol
574 else
575 rtol_ = eps
576 end if
577
578 if (prcsn == 6) then
579 normmax = 1.0e4_dp
580 end if
581 if (present(maxxnorm)) then
582 maxxnorm_ = min(maxxnorm, one/eps)
583 else
584 maxxnorm_ = normmax
585 end if
586
587 if (present(trancond)) then
588 trancond_ = min(trancond, normmax)
589 else
590 trancond_ = normmax
591 end if
592
593 if (present(acondlim)) then
594 acondlim_ = min(acondlim, epsinv)
595 else
596 acondlim_ = epsinv
597 end if
598
599 if (present(msolve)) then
600 precon_ = .true.
601 else
602 precon_ = .false.
603 end if
604
605 !------------------------------------------------------------------
606 ! Print heading and initialize.
607 !------------------------------------------------------------------
608 nprint = min(n,20)
609 !debug = .true.
610 lastiter = .false.
611 flag0 = 0
612 istop_ = flag0
613 beta1 = dnrm2(n, b, 1)
614 ieps = 0.1_dp/eps
615 itn_ = 0
616 qlpiter = 0
617 xnorm_ = zero
618 xl2norm = zero
619 axnorm = zero
620 anorm_ = zero
621 acond_ = one
622 pnorm = zero
623 relresl = zero
624 relaresl = zero
625 x = zero
626 xl2 = zero
627 x1last = x(1)
628
629 if (nout_ > 0) then
630 write(nout_, headerstr) enter, n, beta1, precon_, itnlim_, rtol_, &
631 shift_, maxxnorm_, acondlim_, trancond_
632 end if
633
634 !------------------------------------------------------------------
635 ! Set up y and v for the first Lanczos vector v1.
636 ! y = beta1 P'v1, where P = C**(-1).
637 ! v is really P'v1.
638 !------------------------------------------------------------------
639 y = b
640 r1 = b
641 if ( precon_ ) then
642 call msolve( n, b, y )
643 end if
644
645 beta1 = dot_product(b, y)
646
647 if (beta1 < zero .and. dnrm2(n, y, 1) > eps) then ! M must be indefinite.
648 istop_ = 11
649 end if
650
651 if (beta1 == zero) then ! b = 0 exactly. Stop with x = 0.
652 istop_ = 3
653 end if
654
655 beta1 = sqrt( beta1 ) ! Normalize y to get v1 later.
656
657 if (debug) then
658 write(*,ddebugstr1) ' y_', itn_, ' = ', (y(j), j=1,nprint)
659 write(*,*) 'beta1 ', beta1
660 end if
661
662 !------------------------------------------------------------------
663 ! See if Msolve is symmetric.
664 !------------------------------------------------------------------
665 if (checka_ .and. precon_) then
666 call msolve( n, y, r2 )
667 s = dot_product( y, y )
668 t = dot_product(r1, r2)
669 z = abs( s - t )
670 epsa = (abs(s) + eps) * eps**0.33333_dp
671 if (z > epsa) then
672 istop_ = 10
673 end if
674 end if
675
676 !------------------------------------------------------------------
677 ! See if Aprod is symmetric.
678 !------------------------------------------------------------------
679 if (checka_) then
680 call aprod ( n, y, w ) ! w = A*y
681 call aprod ( n, w, r2 ) ! r2 = A*w
682 s = dot_product( w, w )
683 t = dot_product( y, r2)
684 z = abs( s - t )
685 epsa = (abs(s) + eps) * eps**0.33333_dp
686 if (z > epsa) then
687 istop_ = 9
688 end if
689 end if
690
691 !------------------------------------------------------------------
692 ! Initialize other quantities.
693 !------------------------------------------------------------------
694 tau = zero
695 taul = zero
696 gmin = zero
697 beta = zero
698 betan = beta1
699 dbar = zero
700 dltan = zero
701 eta = zero
702 etal = zero
703 etal2 = zero
704 vepln = zero
705 veplnl = zero
706 veplnl2= zero
707 eplnn = zero
708 gama = zero
709 gamal = zero
710 gamal2 = zero
711 phi = beta1
712 cr1 = -one
713 sr1 = zero
714 cr2 = -one
715 sr2 = zero
716 cs = -one
717 sn = zero
718 ul3 = zero
719 ul2 = zero
720 ul = zero
721 u = zero
722 lines = 1
723 headlines = 30 * lines
724 rnorm_ = betan
725 relres = rnorm_ / (anorm_*xnorm_ + beta1)
726 relares= zero
727 r2 = b
728 w = zero ! vector of zeros
729 wl = zero
730 done = .false.
731
732 if (debug) then
733 write(*,*)
734 write(*,*) 'Checking variable values before main loop'
735 write(*,*) 'istop ', istop_, ' done ', done, ' itn ', itn_, ' QLPiter' , qlpiter
736 write(*,*) 'lastiter', lastiter, ' lines ', lines, ' headlines ', headlines
737 write(*,*) 'beta ', beta, ' tau ', tau, ' taul ', taul, ' phi ', phi
738 write(*,*) 'betan ', betan, ' gmin ', gmin, ' cs ', cs, ' sn ', sn
739 write(*,*) 'cr1 ', cr1, ' sr1 ', sr1, ' cr2 ', cr2, ' sr2 ', sr2
740 write(*,*) 'dltan ', dltan, ' eplnn ', eplnn, ' gama ', gama, ' gamal ', gamal
741 write(*,*) 'gamal2 ', gamal2, ' eta ', eta, ' etal ', etal, 'etal2', etal2
742 write(*,*) 'vepln ', vepln, ' veplnl', veplnl, ' veplnl2 ', veplnl2, ' ul3 ', ul3
743 write(*,*) 'ul2 ', ul2, ' ul ', ul, ' u ', u, ' rnorm ', rnorm_
744 write(*,*) 'xnorm ', xnorm_, ' xl2norm ', xl2norm, ' pnorm ', pnorm, ' Axnorm ', axnorm
745 write(*,*) 'Anorm ', anorm_, ' Acond ', acond_, ' relres ', relres
746 write(*,ddebugstr1) 'w_', itn_-1, ' = ', (wl(j), j=1,nprint)
747 write(*,ddebugstr1) 'w_', itn_, ' = ', (w(j), j=1,nprint)
748 write(*,ddebugstr1) 'x_', itn_, ' = ', (x(j), j=1,nprint)
749 end if
750
751
752 !------------------------------------------------------------------
753 ! Main iteration loop.
754 !------------------------------------------------------------------
755 do while (istop_ <= flag0)
756 itn_ = itn_ + 1 ! k = itn = 1 first time through
757
758 !-----------------------------------------------------------------
759 ! Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
760 ! The general iteration is similar to the case k = 1 with v0 = 0:
761 !
762 ! p1 = Operator * v1 - beta1 * v0,
763 ! alpha1 = v1'p1,
764 ! q2 = p2 - alpha1 * v1,
765 ! beta2^2 = q2'q2,
766 ! v2 = (1/beta2) q2.
767 !
768 ! Again, y = betak P vk, where P = C**(-1).
769 ! .... more description needed.
770 !-----------------------------------------------------------------
771
772 betal = beta; ! betal = betak
773 beta = betan;
774 s = one / beta ! Normalize previous vector (in y).
775 v = s*y; ! v = vk if P = I.
776 call aprod ( n, v, y )
777 if (shift_ /= zero) then
778 y = y - shift_ * v
779 end if
780 if (itn_ >= 2) then
781 y = y + (- beta/betal) * r1
782 end if
783 alfa = dot_product(v, y) ! alphak
784 y = y + (- alfa/beta) * r2
785 r1 = r2
786 r2 = y
787
788 if ( .not. precon_ ) then
789 betan = dnrm2(n, y, 1) ! betan = ||y||_2
790 else
791 call msolve( n, r2, y )
792 betan = dot_product(r2, y) ! betan = betak+1^2
793 if (betan > zero) then
794 betan = sqrt(betan)
795 elseif ( dnrm2(n, y, 1) > eps ) then ! M must be indefinite.
796 istop_ = 11
797 exit
798 end if
799 end if
800
801 if (itn_ == 1) then
802 vec2(1) = alfa
803 vec2(2) = betan
804 pnorm = dnrm2(2, vec2, 1)
805 else
806 vec3(1) = beta
807 vec3(2) = alfa
808 vec3(3) = betan
809 pnorm = dnrm2(3, vec3, 1)
810 end if
811
812
813 if (debug) then
814 write(*,*)
815 write(*,*) 'Lanczos iteration ', itn_
816 write(*,ddebugstr1) 'v_', itn_, ' = ', (v(j), j=1,nprint)
817 write(*,ddebugstr1) 'r1_', itn_, ' = ', (r1(j), j=1,nprint)
818 write(*,ddebugstr1) 'r2_', itn_, ' = ', (r2(j), j=1,nprint)
819 write(*,ddebugstr1) 'y_', itn_, ' = ', (y(j), j=1,nprint)
820 write(*,ddebugstr2) 'alpha_', itn_, ' = ', alfa, ', ', &
821 'beta_', itn_, ' = ', beta, ', ', &
822 'beta_', itn_+1, ' = ', betan, ', ', &
823 'pnorm_', itn_, ' = ', pnorm
824 end if
825
826 ! Apply previous left reflection Qk-1 to get
827 ! [deltak epslnk+1] = [cs sn][dbark 0 ]
828 ! [gbar k dbar k+1] [sn -cs][alfak betak+1].
829
830 dbar = dltan
831 dlta = cs * dbar + sn * alfa ! dlta1 = 0 deltak
832 epln = eplnn;
833 gbar = sn * dbar - cs * alfa ! gbar 1 = alfa1 gbar k
834 eplnn = sn * betan ! eplnn2 = 0 epslnk+1
835 dltan = - cs * betan ! dbar 2 = beta2 dbar k+1
836 dlta_qlp = dlta;
837
838 if (debug) then
839 write(*,*)
840 write(*,*) 'Apply previous left reflection Q_{', itn_-1, ',', itn_, '}'
841 write(*,ddebugstr2) 'c_', itn_-1, ' = ', cs, ', ', &
842 's_', itn_-1, ' = ', sn
843 write(*,ddebugstr2) 'dlta_', itn_, ' = ', dlta, ', ', &
844 'gbar_', itn_, ' = ', gbar, ', ', &
845 'epln_', itn_+1, ' = ', eplnn, ', ', &
846 'dbar_', itn_+1, ' = ', dltan
847 end if
848
849 ! Compute the current left reflection Qk
850 gamal3 = gamal2
851 gamal2 = gamal
852 gamal = gama
853 call symortho(gbar, betan, cs, sn, gama)
854 gama_tmp = gama;
855 taul2 = taul
856 taul = tau
857 tau = cs * phi
858 phi = sn * phi ! phik
859 axnorm = sqrt( axnorm**2 + tau**2 );
860
861 if (debug) then
862 write(*,*)
863 write(*,*) 'Compute the current left reflection Q_{', itn_, ',', itn_+1, '}'
864 write(*,ddebugstr2) 'c_', itn_, ' = ', cs, ', ', &
865 's_', itn_, ' = ', sn
866 write(*,ddebugstr2) 'tau_', itn_, ' = ', tau, ', ', &
867 'phi_', itn_, ' = ', phi, ', ', &
868 'gama_', itn_, ' = ', gama
869 end if
870
871 ! Apply the previous right reflection P{k-2,k}
872
873 if (itn_ > 2) then
874 veplnl2 = veplnl
875 etal2 = etal
876 etal = eta
877 dlta_tmp = sr2 * vepln - cr2 * dlta
878 veplnl = cr2 * vepln + sr2 * dlta
879 dlta = dlta_tmp
880 eta = sr2 * gama
881 gama = -cr2 * gama
882 if (debug) then
883 write(*,*)
884 write(*,*) 'Apply the previous right reflection P_{', itn_-2, ',', itn_, '}'
885 write(*,ddebugstr2) 'cr2_', itn_, ' = ', cr2, ', ', &
886 'sr2_', itn_, ' = ', sr2
887 write(*,ddebugstr2) 'gamal2_', itn_, ' = ', gamal2, ', ', &
888 'gamal_', itn_, ' = ', gamal, ', ', &
889 'gama_', itn_, ' = ', gama
890 write(*,ddebugstr2) 'dlta_', itn_, ' = ', dlta, ', ', &
891 'vepln_', itn_-1, ' = ', veplnl
892 end if
893 end if
894
895
896 ! Compute the current right reflection P{k-1,k}, P_12, P_23,...
897 if (itn_ > 1) then
898 call symortho(gamal, dlta, cr1, sr1, gamal_tmp)
899 gamal = gamal_tmp
900 vepln = sr1 * gama
901 gama = -cr1 * gama
902
903 if (debug) then
904 write(*,*)
905 write(*,*) 'Apply the current right reflection P_{', itn_-1, ',', itn_, '}'
906 write(*,ddebugstr2) 'cr1_ ', itn_, ' = ', cr1, ', ', &
907 'sr1_' , itn_, ' = ', sr1
908 write(*,ddebugstr2) 'gama_', itn_-2, ' = ', gamal2, ', ', &
909 'gama_', itn_-1, ' = ', gamal, ', ', &
910 'gama_', itn_, ' = ', gama
911 write(*,ddebugstr2) 'dlta_', itn_, ' = ', dlta, ', ', &
912 'vepln_', itn_-1, ' = ', veplnl, ', ', &
913 'eta_', itn_, ' = ', eta
914 end if
915 end if
916
917 ! Update xnorm
918
919 xnorml = xnorm_
920 ul4 = ul3
921 ul3 = ul2
922
923 if (itn_ > 2) then
924 ul2 = ( taul2 - etal2 * ul4 - veplnl2 * ul3 ) / gamal2
925 if (debug) then
926 write(*,ddebugstr2) 'tau_', itn_-2, ' = ', taul2, ', ', &
927 'eta_', itn_-2, ' = ', etal2, ', ', &
928 'vepln_', itn_-2, ' = ', veplnl2, ', ', &
929 'gama_', itn_-2, ' = ', gamal2
930 end if
931 end if
932 if (itn_ > 1) then
933 ul = ( taul - etal * ul3 - veplnl * ul2) / gamal
934 if (debug) then
935 write(*,ddebugstr2) 'tau_', itn_-1, ' = ', taul, ', ', &
936 'eta_', itn_-1, ' = ', etal, ', ', &
937 'vepln_', itn_-1, ' = ', veplnl, ', ', &
938 'gamal_', itn_-1, ' = ', gamal
939 end if
940 end if
941
942 vec3(1) = xl2norm
943 vec3(2) = ul2
944 vec3(3) = ul
945 xnorm_tmp = dnrm2(3, vec3, 1) ! norm([xl2norm ul2 ul]);
946 if (abs(gama) > eps) then ! .and. xnorm_tmp < maxxnorm_) then
947 if (debug) then
948 write(*,ddebugstr2) 'tau_', itn_, ' = ', tau, ', ', &
949 'eta_', itn_, ' = ', eta, ', ', &
950 'vepln_', itn_, ' = ', vepln, ', ', &
951 'gama_', itn_, ' = ', gama
952 end if
953 u = (tau - eta*ul2 - vepln*ul) / gama
954 likels = relaresl < relresl
955 vec2(1) = xnorm_tmp
956 vec2(2) = u
957 if (likels .and. dnrm2(2, vec2, 1) > maxxnorm_) then
958 u = zero
959 istop_ = 12
960 end if
961 else
962 u = zero
963 istop_ = 14
964 end if
965 vec2(1) = xl2norm
966 vec2(2) = ul2
967 xl2norm = dnrm2(2, vec2, 1)
968 vec3(1) = xl2norm
969 vec3(2) = ul
970 vec3(3) = u
971 xnorm_ = dnrm2(3, vec3, 1)
972
973 if (acond_ < trancond_ .and. istop_ == flag0 .and. qlpiter == 0) then !! MINRES updates
974 wl2 = wl
975 wl = w
976 if (gama_tmp > eps) then
977 s = one / gama_tmp
978 w = (v - epln*wl2 - dlta_qlp*wl) * s
979 end if
980
981 if (xnorm_ < maxxnorm_) then
982 x1last = x(1)
983 x = x + tau*w
984 else
985 istop_ = 12
986 lastiter = .true.
987 end if
988
989 if (debug) then
990 write(*,*)
991 write(*,*) 'MINRES updates'
992 write(*,ddebugstr2) 'gama_tmp_', itn_, ' = ', gama_tmp, ', ', &
993 'tau_', itn_, ' = ', tau, ', ', &
994 'epln_', itn_, ' = ', epln, ', ', &
995 'dlta_QLP_', itn_, ' = ', dlta_qlp
996 write(*,ddebugstr1) 'v_', itn_ , ' = ', (v(j), j=1,nprint)
997 write(*,ddebugstr1) 'w_', itn_ , ' = ', (w(j), j=1,nprint)
998 end if
999 else !! MINRES-QLP updates
1000 qlpiter = qlpiter + 1;
1001 if (qlpiter == 1) then
1002 xl2 = zero ! vector
1003 if (itn_ > 1) then ! construct w_{k-3}, w_{k-2}, w_{k-1}
1004 if (itn_ > 3) then
1005 wl2 = gamal3*wl2 + veplnl2*wl + etal*w
1006 end if ! w_{k-3}
1007 if (itn_ > 2) then
1008 wl = gamal_qlp*wl + vepln_qlp*w
1009 end if ! w_{k-2}
1010 w = gama_qlp*w
1011 xl2 = x - ul_qlp*wl - u_qlp*w
1012 end if
1013 end if
1014
1015 if (itn_ == 1) then
1016 wl2 = wl
1017 wl = sr1*v
1018 w = -cr1*v
1019 else if (itn_ == 2) then
1020 wl2 = wl
1021 wl = cr1*w + sr1*v
1022 w = sr1*w - cr1*v
1023 else
1024 wl2 = wl
1025 wl = w
1026 w = sr2*wl2 - cr2*v
1027 wl2 = cr2*wl2 + sr2*v
1028 v = cr1*wl + sr1*w
1029 w = sr1*wl - cr1*w
1030 wl = v
1031 end if
1032 x1last = x(1)
1033 xl2 = xl2 + ul2*wl2
1034 x = xl2 + ul *wl + u*w
1035
1036 if (debug) then
1037 write(*,*)
1038 write(*,*) 'MINRESQLP updates'
1039 end if
1040 end if
1041
1042 if (debug) then
1043 write(*,*)
1044 write(*,*) 'Update u, w, x and xnorm'
1045 write(*,ddebugstr2) 'u_', itn_-2, ' = ', ul2, ', ', &
1046 'u_', itn_-1, ' = ', ul, ', ', &
1047 'u_', itn_, ' = ', u
1048 write(*,ddebugstr1) 'w_', itn_-2, ' = ', (wl2(j), j=1,nprint)
1049 write(*,ddebugstr1) 'w_', itn_-1, ' = ', (wl(j), j=1,nprint)
1050 write(*,ddebugstr1) 'w_', itn_, ' = ', (w(j), j=1,nprint)
1051 write(*,ddebugstr1) 'x_', itn_, ' = ', (x(j), j=1,nprint)
1052 write(*,ddebugstr2) 'xnorm_', itn_-2, ' = ', xl2norm, ', ', &
1053 'xnorm_', itn_-1, ' = ', xnorml, ', ', &
1054 'xnorm_', itn_, ' = ', xnorm_
1055 end if
1056
1057 ! Compute the next right reflection P{k-1,k+1}
1058
1059 if (debug) then
1060 write(*,*)
1061 write(*,*) 'Compute the next right reflection P{', itn_-1, itn_+1,'}'
1062 write(*,ddebugstr2) 'gama_', itn_-1, ' = ', gamal, ', ', &
1063 'epln_', itn_+1, ' = ', eplnn
1064 end if
1065
1066 gamal_tmp = gamal
1067 call symortho(gamal_tmp, eplnn, cr2, sr2, gamal)
1068
1069 if (debug) then
1070 write(*,ddebugstr2) 'cr2_', itn_+1, ' = ', cr2, ', ', &
1071 'sr2_', itn_+1, ' = ', sr2, ', ', &
1072 'gama_', itn_-1, ' = ', gamal
1073 end if
1074
1075 ! Store quantities for transfering from MINRES to MINRES-QLP
1076
1077 gamal_qlp = gamal_tmp
1078 vepln_qlp = vepln
1079 gama_qlp = gama
1080 ul2_qlp = ul2
1081 ul_qlp = ul
1082 u_qlp = u
1083
1084 if (debug) then
1085 write(*,*)
1086 write(*,*) 'Store quantities for transfering from MINRES to MINRES-QLP '
1087 write(*,ddebugstr2) 'gama_QLP_', itn_-1, ' = ', gamal_qlp, ', ', &
1088 'vepln_QLP_',itn_, ' = ', vepln_qlp, ', ', &
1089 'gama_QLP_', itn_, ' = ', gama_qlp
1090 write(*,ddebugstr2) 'u_QLP_', itn_-2, ' = ', ul2_qlp, ', ', &
1091 'u_QLP_', itn_-1, ' = ', ul_qlp, ', ', &
1092 'u_QLP_', itn_, ' = ', u_qlp
1093 end if
1094
1095 ! Estimate various norms
1096
1097 abs_gama = abs(gama)
1098 anorml = anorm_
1099 anorm_ = max(anorm_, pnorm, gamal, abs_gama)
1100 if (itn_ == 1) then
1101 gmin = gama
1102 gminl = gmin
1103 else if (itn_ > 1) then
1104 gminl2 = gminl
1105 gminl = gmin
1106 vec3(1) = gminl2
1107 vec3(2) = gamal
1108 vec3(3) = abs_gama
1109 gmin = min(gminl2, gamal, abs_gama)
1110 end if
1111 acondl = acond_
1112 acond_ = anorm_ / gmin
1113 rnorml = rnorm_
1114 relresl = relres
1115 if (istop_ /= 14) rnorm_ = phi
1116 relres = rnorm_ / (anorm_ * xnorm_ + beta1)
1117 vec2(1) = gbar
1118 vec2(2) = dltan
1119 rootl = dnrm2(2, vec2, 1)
1120 arnorml = rnorml * rootl
1121 relaresl = rootl / anorm_
1122
1123 if (debug) then
1124 write(*,*)
1125 write(*,*) 'Estimate various norms '
1126 write(*,ddebugstr2) 'gmin_', itn_, ' = ', gmin, ', ', &
1127 'pnorm_', itn_, ' = ', pnorm, ', ', &
1128 'rnorm_', itn_, ' = ', rnorm_, ', ', &
1129 'Arnorm_', itn_-1, ' = ', arnorml
1130 write(*,ddebugstr2) 'Axnorm_', itn_, ' = ', axnorm, ', ', &
1131 'Anorm_', itn_, ' = ', anorm_, ', ', &
1132 'Acond_', itn_, ' = ', acond_
1133 end if
1134
1135 ! See if any of the stopping criteria are satisfied.
1136
1137 epsx = anorm_*xnorm_*eps
1138 if (istop_ == flag0 .or. istop_ == 14) then
1139 t1 = one + relres
1140 t2 = one + relaresl
1141 end if
1142 if (t1 <= one ) then
1143 istop_ = 5 ! Accurate Ax=b solution
1144 else if (t2 <= one ) then
1145 istop_ = 7 ! Accurate LS solution
1146 else if (relres <= rtol_ ) then
1147 istop_ = 4 ! Good enough Ax=b solution
1148 else if (relaresl <= rtol_ ) then
1149 istop_ = 6 ! Good enough LS solution
1150 else if (epsx >= beta1 ) then
1151 istop_ = 2 ! x is an eigenvector
1152 else if (xnorm_ >= maxxnorm_) then
1153 istop_ = 12 ! xnorm exceeded its limit
1154 else if (acond_ >= acondlim_ .or. acond_ >= ieps) then
1155 istop_ = 13 ! Huge Acond
1156 else if (itn_ >= itnlim_ ) then
1157 istop_ = 8 ! Too many itns
1158 else if (betan < eps ) then
1159 istop_ = 1 ! Last iteration of Lanczos, rarely happens
1160 end if
1161
1162 if (disable_ .and. itn_ < itnlim_) then
1163 istop_ = flag0
1164 done = .false.
1165 if (axnorm < rtol_*anorm_*xnorm_) then
1166 istop_ = 15
1167 lastiter = .false.
1168 end if
1169 end if
1170
1171 if (istop_ /= flag0) then
1172 done = .true.
1173 if (istop_ == 6 .or. istop_ == 7 .or. istop_ == 12 .or. istop_ == 13) then
1174 lastiter = .true.
1175 end if
1176 if (lastiter) then
1177 itn_ = itn_ - 1
1178 acond_ = acondl
1179 rnorm_ = rnorml
1180 relres = relresl
1181 end if
1182
1183 call aprod ( n, x, r1 )
1184 r1 = b - r1 + shift_*x ! r1 to temporarily store residual vector
1185 call aprod ( n, r1, wl2 ) ! wl2 to temporarily store A*r1
1186 wl2 = wl2 - shift_*r1
1187 arnorm_ = dnrm2(n, wl2, 1)
1188 if (rnorm_ > zero .and. anorm_ > zero) then
1189 relares = arnorm_ / (anorm_*rnorm_)
1190 end if
1191 end if
1192
1193 if (nout_ > 0 .and. .not. lastiter .and. mod(itn_-1,lines) == 0) then
1194 if (itn_ == 101) then
1195 lines = 10
1196 headlines = 30*lines
1197 else if (itn_ == 1001) then
1198 lines = 100
1199 headlines = 30*lines
1200 end if
1201
1202 if (qlpiter == 1) then
1203 qlpstr = ' P'
1204 else
1205 qlpstr = ' '
1206 end if
1207 end if
1208
1209
1210 ! See if it is time to print something.
1211
1212 if (nout_ > 0) then
1213 prnt = .false.
1214 if (n <= 40 ) prnt = .true.
1215 if (itn_ <= 10 ) prnt = .true.
1216 if (itn_ >= itnlim_ - 10) prnt = .true.
1217 if (mod(itn_-1,10) == 0 ) prnt = .true.
1218 if (qlpiter == 1 ) prnt = .true.
1219 if (acond_ >= 0.01_dp/eps ) prnt = .true.
1220 if (istop_ /= flag0 ) prnt = .true.
1221
1222 if ( prnt ) then
1223 if (itn_ == 1) write(nout_, tableheaderstr)
1224 write(nout_, itnstr) itn_-1, x1last, xnorml, &
1225 rnorml, arnorml, relresl, relaresl, anorml, acondl, qlpstr
1226 if (mod(itn_,10) == 0) write(nout_, "(1x)")
1227 end if
1228 end if
1229
1230 if (debug) then
1231 write(*,*) 'istop = ', istop_
1232 end if
1233 if (istop_ /= flag0) exit
1234 enddo
1235 !===================================================================
1236 ! End of iteration loop.
1237 !===================================================================
1238
1239 ! Optional outputs
1240
1241 if (present(istop)) then
1242 istop = istop_
1243 end if
1244
1245 if (present(itn)) then
1246 itn = itn_
1247 end if
1248
1249 if (present(rnorm)) then
1250 rnorm = rnorm_
1251 end if
1252
1253 if (present(arnorm)) then
1254 arnorm = arnorm_
1255 end if
1256
1257 if (present(xnorm)) then
1258 xnorm = xnorm_
1259 end if
1260
1261 if (present(anorm)) then
1262 anorm = anorm_
1263 end if
1264
1265 if (present(acond)) then
1266 acond = acond_
1267 end if
1268
1269 if ( prnt ) then
1270 write(nout_, itnstr) itn_, x(1), xnorm_, &
1271 rnorm_, arnorm_, relres, relares, anorm_, acond_, qlpstr
1272 end if
1273
1274 ! Display final status.
1275
1276 if (nout_ > 0) then
1277 write(nout_, finalstr1) exitt, 'istop =', istop_, 'itn =', itn_, &
1278 exitt, 'Anorm =', anorm_, 'Acond =', acond_, &
1279 exitt, 'rnorm =', rnorm_, 'Arnorm =', arnorm_, &
1280 exitt, 'xnorm =', xnorm_
1281 write(nout_, finalstr2) exitt, msg(istop_)
1282 end if
1283
1284 return
1285end subroutine minresqlp
1286
1287 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1288
1289 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1332 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1333
1334 subroutine symortho(a, b, c, s, r)
1335
1336 real(dp), intent(in) :: a, b
1337 real(dp), intent(out) :: c, s, r
1338
1339 !obsolete intrinsic :: abs, sqrt
1340
1341 ! Local variables
1342 logical, parameter :: debug = .false.
1343 real(dp) :: t
1344 real(dp) :: abs_a, abs_b
1345
1346 abs_a = abs(a)
1347 abs_b = abs(b)
1348
1349 if (abs_b <= realmin) then
1350 s = zero
1351 r = abs_a
1352 if (a == zero) then
1353 c = one
1354 else
1355 c = a / abs_a
1356 end if
1357
1358 else if (abs_a <= realmin) then
1359 c = zero;
1360 r = abs_b
1361 s = b / abs_b
1362
1363 else if (abs_b > abs_a) then
1364 t = a / b
1365 s = (b / abs_b) / sqrt(one + t**2)
1366 c = s * t
1367 r = b / s ! computationally better than r = a / c since |c| <= |s|
1368
1369 else
1370 t = b / a
1371 c = (a / abs_a) / sqrt(one + t**2)
1372 s = c * t;
1373 r = a / c ! computationally better than r = b / s since |s| <= |c|
1374 end if
1375
1376 if (debug) then
1377 write(*,*) 'c = ', c, ', s = ', s, ', r = ', r
1378 end if
1379
1380 end subroutine symortho
1381
1382end module minresqlpmodule
real(dp) function, public dnrm2(n, x, incx)
DNRM2 returns the euclidean norm of a vector.
Defines precision and range in real(kind=dp) and integer(kind=ip) for portability and a few constants...
real(dp), parameter, public eps
real(dp), parameter, public one
integer, parameter, public prcsn
real(dp), parameter, public zero
integer, parameter, public ip
integer, parameter, public dp
real(dp), parameter, public realmin
MINRESQLP solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite a...
subroutine, public minresqlp(n, Aprod, b, shift, Msolve, disable, nout, itnlim, rtol, maxxnorm, trancond, Acondlim, x, istop, itn, rnorm, Arnorm, xnorm, Anorm, Acond)
Solution of linear equation system or least squares problem.
subroutine, public symortho(a, b, c, s, r)
SymOrtho: Stable Householder reflection.