37 REAL(
mpd),
DIMENSION(:),
ALLOCATABLE ::
matv
39 REAL(
mpd),
DIMENSION(:),
ALLOCATABLE ::
matl
40 REAL(
mpd),
DIMENSION(:),
ALLOCATABLE ::
vecn
41 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE ::
nparblock
42 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE ::
ioffblock
43 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE ::
ioffrow
44 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE ::
ioffpar
62 INTEGER(mpl) :: length
64 INTEGER(mpi),
INTENT(IN) :: n
65 INTEGER(mpi),
INTENT(IN) :: m
66 INTEGER(mpi),
INTENT(IN) :: l
67 INTEGER(mpl),
INTENT(IN) :: s
68 INTEGER(mpi),
INTENT(IN) :: k
129 INTEGER(mpl) :: ioff1
130 INTEGER(mpl) :: ioff2
131 INTEGER(mpl) :: ioff3
134 INTEGER(mpl) :: length
138 REAL(mpd),
INTENT(IN) :: a(matsize)
159 ioff1=int(k-1,mpl)*int(
npar,mpl)
162 nrm = sqrt(dot_product(
vecn(1:kn),
vecn(1:kn)))
163 IF (nrm == 0.0_mpd) cycle
165 IF (
vecn(kn) >= 0.0_mpd)
THEN
171 nrm = sqrt(dot_product(
vecn(1:kn),
vecn(1:kn)))
176 sp=dot_product(
vecn(1:kn),
matv(ioff2+1:ioff2+kn))
177 matv(ioff2+1:ioff2+kn)=
matv(ioff2+1:ioff2+kn)-2.0_mpd*
vecn(1:kn)*sp
181 ioff3=int(k-1,mpl)*int(
ncon,mpl)
184 matv(ioff1+1:ioff1+kn-1)=
vecn(1:kn-1)
223 INTEGER(mpi) :: ibcon
224 INTEGER(mpi) :: iblast
225 INTEGER(mpi) :: iboff
226 INTEGER(mpi) :: ibpar
227 INTEGER(mpi) :: ifirst
228 INTEGER(mpi) :: ilast
230 INTEGER(mpl) :: ioff1
231 INTEGER(mpl) :: ioff2
232 INTEGER(mpl) :: ioff3
233 INTEGER(mpi) :: iclast
234 INTEGER(mpi) :: icoff
235 INTEGER(mpi) :: iplast
236 INTEGER(mpi) :: ipoff
246 REAL(mpd),
INTENT(IN) :: a(matsize)
247 INTEGER(mpi),
INTENT(IN) :: bpar(2,nblock+1)
248 INTEGER(mpi),
INTENT(IN) :: bcon(3,ncon+1)
249 INTEGER(mpi),
INTENT(IN) :: rcon(4,ncon)
261 DO ibcon=bpar(2,ibpar)+1, bpar(2,ibpar+1)
262 ncb=bcon(1,ibcon+1)-bcon(1,ibcon)
263 npb=bcon(3,ibcon)+1-bcon(2,ibcon)
266 DO i=bcon(1,ibcon),bcon(1,ibcon+1)-1
277 iblast=bpar(1,ibpar+1)
281 npb=int(
ioffrow(k+1)-ioff1,mpi)
285 ioff3=int(k-1,mpl)*int(ncon,mpl)
286 matl(ioff3+iclast-ncol-icoff:ioff3+iclast-icoff)=
matv(ioff1+npb-ncol:ioff1+npb)
290 nparblock(ibpar)=bpar(1,ibpar+1)-bpar(1,ibpar)
296 iblast=bpar(1,ibpar+1)
299 IF(iclast <= icoff) cycle
300 ibcon=bpar(2,ibpar+1)
303 DO k=iclast,icoff+1,-1
316 IF (ifirst > kn) cycle
320 npb=int(
ioffrow(k+1)-ioff1,mpi)
330 vecn(ifirst:ilast)=
matv(ioff1+1:ioff1+ncol)
331 nrm = sqrt(dot_product(
vecn(ifirst:ilast),
vecn(ifirst:ilast)))
332 IF (nrm == 0.0_mpd) cycle
334 IF (
vecn(kn) >= 0.0_mpd)
THEN
341 nrm = sqrt(dot_product(
vecn(ifirst:ilast),
vecn(ifirst:ilast))+
vecn(kn)*
vecn(kn))
342 vecn(ifirst:ilast)=
vecn(ifirst:ilast)/nrm
345 nrm = sqrt(dot_product(
vecn(ifirst:ilast),
vecn(ifirst:ilast)))
346 vecn(ifirst:ilast)=
vecn(ifirst:ilast)/nrm
349 ioff3=int(k1-2,mpl)*int(ncon,mpl)
355 sp=dot_product(
vecn(ifirst:ilast),
matv(ioff2+1:ioff2+ncol))
356 IF (sp /= 0.0_mpd)
THEN
358 matv(ioff2+1:ioff2+ncol)=
matv(ioff2+1:ioff2+ncol)-2.0_mpd*
vecn(ifirst:ilast)*sp
361 matl(ioff3+i-icoff:ioff3+k-icoff)=
matl(ioff3+i-icoff:ioff3+k-icoff)-2.0_mpd*
vecn(in:kn)*sp
373 matv(ioff1+1:ioff1+ncol)=
vecn(ifirst:ilast)
374 matv(ioff1+ncol+1:ioff1+npb)=0.0_mpd
401 INTEGER(mpi) :: icoff
402 INTEGER(mpi) :: iclast
403 INTEGER(mpi) :: ifirst
404 INTEGER(mpi) :: ilast
405 INTEGER(mpl) :: ioff2
410 INTEGER(mpi) :: nconb
411 INTEGER(mpi) :: nparb
414 INTEGER(mpi),
INTENT(IN) :: m
415 REAL(mpd),
INTENT(IN OUT) :: x(INT(npar,mpl)*INT(m,mpl))
416 LOGICAL,
INTENT(IN) :: t
437 sp=dot_product(
vecn(ifirst:ilast),x(ioff2+ifirst:ioff2+ilast))+
vecn(kn)*x(ioff2+kn)
438 x(ioff2+ifirst:ioff2+ilast)=x(ioff2+ifirst:ioff2+ilast)-2.0_mpd*
vecn(ifirst:ilast)*sp
439 x(ioff2+kn)=x(ioff2+kn)-2.0_mpd*
vecn(kn)*sp
463 INTEGER(mpi) :: ifirst
464 INTEGER(mpi) :: ilast
470 INTEGER(mpi),
INTENT(IN) :: m
471 REAL(mpd),
INTENT(IN OUT) :: x(INT(m,mpl)*INT(npar,mpl))
472 LOGICAL,
INTENT(IN) :: t
476 IF (.not.t) k=
ncon+1-j
488 sp=dot_product(
vecn(1:kn),x(i:iend:m))
489 x(i:iend:m)=x(i:iend:m)-2.0_mpd*
vecn(1:kn)*sp
510 INTEGER(mpl) :: ioff2
512 INTEGER(mpi) :: ifirst
513 INTEGER(mpi) :: ilast
519 REAL(mpd),
INTENT(IN OUT) :: x(INT(npar,mpl)*INT(npar,mpl))
520 LOGICAL,
INTENT(IN) :: t
536 iend=int(npar,mpl)*int(kn,mpl)
538 sp=dot_product(
vecn(1:kn),x(i:iend:npar))
539 x(i:iend:npar)=x(i:iend:npar)-2.0_mpd*
vecn(1:kn)*sp
543 sp=dot_product(
vecn(1:kn),x(ioff2+1:ioff2+kn))
544 x(ioff2+1:ioff2+kn)=x(ioff2+1:ioff2+kn)-2.0_mpd*
vecn(1:kn)*sp
571 INTEGER(mpi) :: ibpar
572 INTEGER(mpi) :: icoff
573 INTEGER(mpi) :: iclast
574 INTEGER(mpi) :: ifirst
575 INTEGER(mpi) :: ilast
576 INTEGER(mpi) :: ilasti
577 INTEGER(mpl) :: ioff2
578 INTEGER(mpi) :: ioffp
583 INTEGER(mpl) :: length
584 INTEGER(mpi) :: nconb
585 INTEGER(mpi) :: nparb
587 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: Av
589 INTEGER(mpl),
INTENT(IN) :: s
590 REAL(mpd),
INTENT(IN OUT) :: A(s)
591 INTEGER(mpl),
INTENT(IN) :: roff(npar)
592 LOGICAL,
INTENT(IN) :: t
595 SUBROUTINE aprod(n,l,x,is,ie,y)
597 INTEGER(mpi),
INTENT(in) :: n
598 INTEGER(mpl),
INTENT(in) :: l
599 REAL(mpd),
INTENT(IN) :: x(n)
600 INTEGER(mpi),
INTENT(in) :: is
601 INTEGER(mpi),
INTENT(in) :: ie
602 REAL(mpd),
INTENT(OUT) :: y(n)
608 CALL mpalloc(av,length,
'qlssq: A*v')
632 CALL aprod(nparb,int(ioffp,mpl),
vecn(1:nparb),ifirst,ilast,av(1:nparb))
633 CALL aprod(nparb,int(ioffp,mpl),
vecn(1:nparb),kn,kn,av(1:nparb))
637 vtav=dot_product(
vecn(ifirst:ilast),av(ifirst:ilast))+
vecn(kn)*av(kn)
645 ioff2=roff(i+ioffp)+ioffp
648 a(ioff2+ifirst:ioff2+ilasti)=a(ioff2+ifirst:ioff2+ilasti)+2.0_mpd* &
649 ((2.0_mpd*
vecn(i)*vtav-av(i))*
vecn(ifirst:ilasti))
659 ioff2=roff(i+ioffp)+ioffp
661 a(ioff2+1:ioff2+i)=a(ioff2+1:ioff2+i)-2.0_mpd*av(1:i)*
vecn(i)
665 ioff2=roff(kn+ioffp)+ioffp
666 a(ioff2+kn)=a(ioff2+kn)+2.0_mpd*((2.0_mpd*
vecvk(l)*vtav-av(kn))*
vecvk(l)-av(kn)*
vecvk(l))
669 ioff2=roff(i+ioffp)+ioffp
671 a(ioff2+ifirst:ioff2+ilast)=a(ioff2+ifirst:ioff2+ilast)-2.0_mpd*
vecn(ifirst:ilast)*av(i)
672 a(ioff2+kn)=a(ioff2+kn)-2.0_mpd*
vecvk(l)*av(i)
703 INTEGER(mpi) :: ifirst
704 INTEGER(mpi) :: ilast
705 INTEGER(mpi) :: ifirst2
706 INTEGER(mpi) :: ilast2
707 INTEGER(mpl) :: ioff1
708 INTEGER(mpl) :: ioff2
709 INTEGER(mpi) :: istat(3)
719 INTEGER(mpl) :: length
725 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecAv
726 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matvtvp
727 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matvtAvp
728 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matCoeff
729 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: irangeCoeff
731 INTEGER(mpi),
INTENT(IN) :: m
732 REAL(mpd),
INTENT(IN OUT) :: B(npar*m-(m*m-m)/2)
733 LOGICAL,
INTENT(IN) :: t
736 SUBROUTINE aprod(n,l,x,is,ie,y)
738 INTEGER(mpi),
INTENT(in) :: n
739 INTEGER(mpl),
INTENT(in) :: l
740 REAL(mpd),
INTENT(IN) :: x(n)
741 INTEGER(mpi),
INTENT(in) :: is
742 INTEGER(mpi),
INTENT(in) :: ie
743 REAL(mpd),
INTENT(OUT) :: y(n)
749 CALL mpalloc(vecav,length,
'qlpssq: A*v')
751 CALL mpalloc(matvtvp,length,
"qlpssq: v^t*v'")
753 CALL mpalloc(matvtavp,length,
"qlpssq: v^t*(A*v')")
755 CALL mpalloc(matcoeff,length,
'qlpssq: coefficients')
758 CALL mpalloc(irangecoeff,2_mpl,length,
'qlpssq: non zero coefficient range')
768 ioff1=int(k-1,mpl)*int(
ncon,mpl)
769 irangecoeff(1,k)=
ncon
780 CALL aprod(npar,0_mpl,
vecn(1:npar),ifirst,ilast,vecav(1:npar))
781 CALL aprod(npar,0_mpl,
vecn(1:npar),kn,kn,vecav(1:npar))
787 ioff2=int(k2-1,mpl)*int(
ncon,mpl)
792 IF (kn >= ifirst2.AND.kn <= ilast2) v2kn=
matv(
ioffrow(k2)+1+kn-ifirst2)
794 l1=max(ifirst,ifirst2)
797 IF (l1 <= l2) vtvp=vtvp+dot_product(
vecn(l1:l2), &
800 IF (abs(vtvp) > 16.0_mpd*epsilon(vtvp))
THEN
801 matvtvp(ioff1+k2)=vtvp
802 matvtvp(ioff2+k)=vtvp
805 matvtvp(ioff1+k)=1.0_mpd
815 matvtavp(ioff1+k2)=
vecvk(k2)*vecav(kn2)+dot_product(vecav(ifirst2:ilast2), &
820 vtav=matvtavp(ioff1+k)
825 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*
vecn(i)*vtav-vecav(i))*
vecn(l)-vecav(l)*
vecn(i))
833 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*
vecn(l)
843 ioff1=int(k-1,mpl)*int(
ncon,mpl)
856 vtav=matvtavp(ioff1+k)+dot_product(matcoeff(ioff1+l1:ioff1+l2),matvtvp(ioff1+l1:ioff1+l2))
861 IF (matcoeff(ioff1+k2) == 0._mpd) cycle
862 if (istat(1)==0) istat(1)=k2
870 vecav(ifirst2:ilast2)=vecav(ifirst2:ilast2)+matcoeff(ioff1+k2)* &
872 vecav(kn2)=vecav(kn2)+matcoeff(ioff1+k2)*
vecvk(k2)
880 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*
vecn(i)*vtav-vecav(i))*
vecn(l)-vecav(l)*
vecn(i))
888 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*
vecn(l)
896 ioff2=int(k2-1,mpl)*int(
ncon,mpl)
897 vtvp=matvtvp(ioff1+k2)
901 vtavp=matvtavp(ioff2+k)
902 IF (l1 <= l2) vtavp=vtavp+dot_product(matcoeff(ioff2+l1:ioff2+l2),matvtvp(ioff1+l1:ioff1+l2))
905 matcoeff(ioff2+k)=matcoeff(ioff2+k)+2.0_mpd*(2.0_mpd*vtav*vtvp-vtavp)
906 IF (vtvp /= 0._mpd)
THEN
907 l1=min(l1,irangecoeff(1,k))
908 l2=max(l2,irangecoeff(2,k))
909 matcoeff(ioff2+l1:ioff2+l2)=matcoeff(ioff2+l1:ioff2+l2)-2.0_mpd*matcoeff(ioff1+l1:ioff1+l2)*vtvp
938 INTEGER(mpi) :: ibpar
939 INTEGER(mpi) :: icoff
940 INTEGER(mpi) :: iclast
941 INTEGER(mpl) :: idiag
943 REAL(mpd),
INTENT(OUT) :: emin
944 REAL(mpd),
INTENT(OUT) :: emax
951 idiag=int(
ncon,mpl)*int(icoff,mpl)+1
953 IF (abs(emax) < abs(
matl(idiag))) emax=
matl(idiag)
954 IF (abs(emin) > abs(
matl(idiag))) emin=
matl(idiag)
973 INTEGER(mpi) :: icoff
974 INTEGER(mpi) :: iclast
975 INTEGER(mpl) :: idiag
977 INTEGER(mpi) :: nconb
979 REAL(mpd),
INTENT(IN) :: d(ncon)
980 REAL(mpd),
INTENT(OUT) :: y(ncon)
986 idiag=int(ncon,mpl)*int(iclast-1,mpl)+nconb
988 y(k)=(d(k)-dot_product(
matl(idiag+1:idiag+nconb-k),y(k+1:nconb)))/
matl(idiag)
1000 INTEGER(mpi),
INTENT(IN) :: ib
1013 INTEGER(mpi) :: ifirst
1014 INTEGER(mpi) :: ilast
1015 INTEGER(mpl) :: ioff1
1016 INTEGER(mpl) :: ioff2
1017 INTEGER(mpi) :: istat(6)
1032 v1=0.;v2=0.;v3=0.;v4=0.
1040 IF (
vecn(j) /= 0.0_mpd)
THEN
1042 IF (istat(3) == 0)
THEN
1052 IF (
matl(ioff2+j) /= 0.0_mpd)
THEN
1054 IF (istat(6) == 0)
THEN
1066100
FORMAT(
" qldump",7i8,5g13.5,2i8)
subroutine monpgs(i)
Progress monitoring.
subroutine qlmrq(x, m, t)
Multiply right by Q(t).
subroutine qldump()
Print statistics.
subroutine qlsmq(x, t)
Similarity transformation by Q(t).
subroutine qlpssq(aprod, B, m, t)
Partial similarity transformation by Q(t).
subroutine qldecb(a, bpar, bcon, rcon)
QL decomposition (for disjoint block matrix).
subroutine qldec(a)
QL decomposition (as single block).
subroutine qlmlq(x, m, t)
Multiply left by Q(t) (per block).
subroutine qlsetb(ib)
Set block.
subroutine qlbsub(d, y)
Backward substitution (per block).
subroutine qlini(n, m, l, s, k)
Initialize QL decomposition.
subroutine qlgete(emin, emax)
Get eigenvalues.
subroutine qlssq(aprod, A, s, roff, t)
Similarity transformation by Q(t).
(De)Allocate vectors and arrays.
integer, parameter mpd
double precision
integer(mpi) ncon
number of constraints
integer(mpi) iblock
active block
real(mpd), dimension(:), allocatable vecvk
secondary diagonal of matV (including last element)
integer(mpi), dimension(:,:), allocatable irangeparnz
range for non zero part (except vecVk)
integer(mpi) monpg
flag for progress monitoring
integer(mpl), dimension(:), allocatable ioffrow
row offsets in matV (for constrint block)
integer(mpl) matsize
size of contraints matrix
real(mpd), dimension(:), allocatable matv
unit normals (v_i) of Householder reflectors
real(mpd), dimension(:), allocatable matl
lower diagonal matrix L
integer(mpi), dimension(:), allocatable ioffpar
parameter number offsets for matV ( " )
integer(mpi) nblock
number of blocks
real(mpd), dimension(:), allocatable vecn
normal vector
integer(mpi) npar
number of parameters
integer(mpi), dimension(:), allocatable ioffblock
block offset (1.
integer(mpi), dimension(:), allocatable nparblock
number of parameters in block