26#ifdef SCOREP_USER_ENABLE
27#include "scorep/SCOREP_User.inc"
41 REAL(
mpd),
DIMENSION(:),
ALLOCATABLE ::
matv
43 REAL(
mpd),
DIMENSION(:),
ALLOCATABLE ::
matl
44 REAL(
mpd),
DIMENSION(:),
ALLOCATABLE ::
vecn
45 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE ::
nparblock
46 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE ::
ioffblock
47 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE ::
ioffrow
48 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE ::
ioffpar
66 INTEGER(mpl) :: length
68 INTEGER(mpi),
INTENT(IN) :: n
69 INTEGER(mpi),
INTENT(IN) :: m
70 INTEGER(mpi),
INTENT(IN) :: l
71 INTEGER(mpl),
INTENT(IN) :: s
72 INTEGER(mpi),
INTENT(IN) :: k
133 INTEGER(mpl) :: ioff1
134 INTEGER(mpl) :: ioff2
135 INTEGER(mpl) :: ioff3
138 INTEGER(mpl) :: length
142 REAL(mpd),
INTENT(IN) :: a(matsize)
163 ioff1=int(k-1,mpl)*int(
npar,mpl)
166 nrm = sqrt(dot_product(
vecn(1:kn),
vecn(1:kn)))
167 IF (nrm == 0.0_mpd) cycle
169 IF (
vecn(kn) >= 0.0_mpd)
THEN
175 nrm = sqrt(dot_product(
vecn(1:kn),
vecn(1:kn)))
180 sp=dot_product(
vecn(1:kn),
matv(ioff2+1:ioff2+kn))
181 matv(ioff2+1:ioff2+kn)=
matv(ioff2+1:ioff2+kn)-2.0_mpd*
vecn(1:kn)*sp
185 ioff3=int(k-1,mpl)*int(
ncon,mpl)
188 matv(ioff1+1:ioff1+kn-1)=
vecn(1:kn-1)
227 INTEGER(mpi) :: ibcon
228 INTEGER(mpi) :: iblast
229 INTEGER(mpi) :: iboff
230 INTEGER(mpi) :: ibpar
231 INTEGER(mpi) :: ifirst
232 INTEGER(mpi) :: ilast
234 INTEGER(mpl) :: ioff1
235 INTEGER(mpl) :: ioff2
236 INTEGER(mpl) :: ioff3
237 INTEGER(mpi) :: iclast
238 INTEGER(mpi) :: icoff
239 INTEGER(mpi) :: iplast
240 INTEGER(mpi) :: ipoff
250 REAL(mpd),
INTENT(IN) :: a(matsize)
251 INTEGER(mpi),
INTENT(IN) :: bpar(2,nblock+1)
252 INTEGER(mpi),
INTENT(IN) :: bcon(3,ncon+1)
253 INTEGER(mpi),
INTENT(IN) :: rcon(4,ncon)
256#ifdef SCOREP_USER_ENABLE
257 scorep_user_region_by_name_begin(
"UR_qldecb", scorep_user_region_type_common)
268 DO ibcon=bpar(2,ibpar)+1, bpar(2,ibpar+1)
269 ncb=bcon(1,ibcon+1)-bcon(1,ibcon)
270 npb=bcon(3,ibcon)+1-bcon(2,ibcon)
273 DO i=bcon(1,ibcon),bcon(1,ibcon+1)-1
284 iblast=bpar(1,ibpar+1)
288 npb=int(
ioffrow(k+1)-ioff1,mpi)
292 ioff3=int(k-1,mpl)*int(ncon,mpl)
293 matl(ioff3+iclast-ncol-icoff:ioff3+iclast-icoff)=
matv(ioff1+npb-ncol:ioff1+npb)
297 nparblock(ibpar)=bpar(1,ibpar+1)-bpar(1,ibpar)
303 iblast=bpar(1,ibpar+1)
306 IF(iclast <= icoff) cycle
307 ibcon=bpar(2,ibpar+1)
310 DO k=iclast,icoff+1,-1
323 IF (ifirst > kn) cycle
327 npb=int(
ioffrow(k+1)-ioff1,mpi)
337 vecn(ifirst:ilast)=
matv(ioff1+1:ioff1+ncol)
338 nrm = sqrt(dot_product(
vecn(ifirst:ilast),
vecn(ifirst:ilast)))
339 IF (nrm == 0.0_mpd) cycle
341 IF (
vecn(kn) >= 0.0_mpd)
THEN
348 nrm = sqrt(dot_product(
vecn(ifirst:ilast),
vecn(ifirst:ilast))+
vecn(kn)*
vecn(kn))
349 vecn(ifirst:ilast)=
vecn(ifirst:ilast)/nrm
352 nrm = sqrt(dot_product(
vecn(ifirst:ilast),
vecn(ifirst:ilast)))
353 vecn(ifirst:ilast)=
vecn(ifirst:ilast)/nrm
356 ioff3=int(k1-2,mpl)*int(ncon,mpl)
362 sp=dot_product(
vecn(ifirst:ilast),
matv(ioff2+1:ioff2+ncol))
363 IF (sp /= 0.0_mpd)
THEN
365 matv(ioff2+1:ioff2+ncol)=
matv(ioff2+1:ioff2+ncol)-2.0_mpd*
vecn(ifirst:ilast)*sp
368 matl(ioff3+i-icoff:ioff3+k-icoff)=
matl(ioff3+i-icoff:ioff3+k-icoff)-2.0_mpd*
vecn(in:kn)*sp
380 matv(ioff1+1:ioff1+ncol)=
vecn(ifirst:ilast)
381 matv(ioff1+ncol+1:ioff1+npb)=0.0_mpd
387#ifdef SCOREP_USER_ENABLE
388 scorep_user_region_by_name_end(
"UR_qldecb")
411 INTEGER(mpi) :: icoff
412 INTEGER(mpi) :: iclast
413 INTEGER(mpi) :: ifirst
414 INTEGER(mpi) :: ilast
415 INTEGER(mpl) :: ioff2
420 INTEGER(mpi) :: nconb
421 INTEGER(mpi) :: nparb
424 INTEGER(mpi),
INTENT(IN) :: m
425 REAL(mpd),
INTENT(IN OUT) :: x(INT(npar,mpl)*INT(m,mpl))
426 LOGICAL,
INTENT(IN) :: t
447 sp=dot_product(
vecn(ifirst:ilast),x(ioff2+ifirst:ioff2+ilast))+
vecn(kn)*x(ioff2+kn)
448 x(ioff2+ifirst:ioff2+ilast)=x(ioff2+ifirst:ioff2+ilast)-2.0_mpd*
vecn(ifirst:ilast)*sp
449 x(ioff2+kn)=x(ioff2+kn)-2.0_mpd*
vecn(kn)*sp
473 INTEGER(mpi) :: ifirst
474 INTEGER(mpi) :: ilast
480 INTEGER(mpi),
INTENT(IN) :: m
481 REAL(mpd),
INTENT(IN OUT) :: x(INT(m,mpl)*INT(npar,mpl))
482 LOGICAL,
INTENT(IN) :: t
486 IF (.not.t) k=
ncon+1-j
498 sp=dot_product(
vecn(1:kn),x(i:iend:m))
499 x(i:iend:m)=x(i:iend:m)-2.0_mpd*
vecn(1:kn)*sp
520 INTEGER(mpl) :: ioff2
522 INTEGER(mpi) :: ifirst
523 INTEGER(mpi) :: ilast
529 REAL(mpd),
INTENT(IN OUT) :: x(INT(npar,mpl)*INT(npar,mpl))
530 LOGICAL,
INTENT(IN) :: t
546 iend=int(npar,mpl)*int(kn,mpl)
548 sp=dot_product(
vecn(1:kn),x(i:iend:npar))
549 x(i:iend:npar)=x(i:iend:npar)-2.0_mpd*
vecn(1:kn)*sp
553 sp=dot_product(
vecn(1:kn),x(ioff2+1:ioff2+kn))
554 x(ioff2+1:ioff2+kn)=x(ioff2+1:ioff2+kn)-2.0_mpd*
vecn(1:kn)*sp
581 INTEGER(mpi) :: ibpar
582 INTEGER(mpi) :: icoff
583 INTEGER(mpi) :: iclast
584 INTEGER(mpi) :: ifirst
585 INTEGER(mpi) :: ilast
586 INTEGER(mpi) :: ilasti
587 INTEGER(mpl) :: ioff2
588 INTEGER(mpi) :: ioffp
593 INTEGER(mpl) :: length
594 INTEGER(mpi) :: nconb
595 INTEGER(mpi) :: nparb
597 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: Av
599 INTEGER(mpl),
INTENT(IN) :: s
600 REAL(mpd),
INTENT(IN OUT) :: A(s)
601 INTEGER(mpl),
INTENT(IN) :: roff(npar)
602 LOGICAL,
INTENT(IN) :: t
605 SUBROUTINE aprod(n,l,x,is,ie,y)
607 INTEGER(mpi),
INTENT(in) :: n
608 INTEGER(mpl),
INTENT(in) :: l
609 REAL(mpd),
INTENT(IN) :: x(n)
610 INTEGER(mpi),
INTENT(in) :: is
611 INTEGER(mpi),
INTENT(in) :: ie
612 REAL(mpd),
INTENT(OUT) :: y(n)
616#ifdef SCOREP_USER_ENABLE
617 scorep_user_region_by_name_begin(
"UR_qlssq", scorep_user_region_type_common)
621 CALL mpalloc(av,length,
'qlssq: A*v')
645 CALL aprod(nparb,int(ioffp,mpl),
vecn(1:nparb),ifirst,ilast,av(1:nparb))
646 CALL aprod(nparb,int(ioffp,mpl),
vecn(1:nparb),kn,kn,av(1:nparb))
650 vtav=dot_product(
vecn(ifirst:ilast),av(ifirst:ilast))+
vecn(kn)*av(kn)
658 ioff2=roff(i+ioffp)+ioffp
661 a(ioff2+ifirst:ioff2+ilasti)=a(ioff2+ifirst:ioff2+ilasti)+2.0_mpd* &
662 ((2.0_mpd*
vecn(i)*vtav-av(i))*
vecn(ifirst:ilasti))
672 ioff2=roff(i+ioffp)+ioffp
674 a(ioff2+1:ioff2+i)=a(ioff2+1:ioff2+i)-2.0_mpd*av(1:i)*
vecn(i)
678 ioff2=roff(kn+ioffp)+ioffp
679 a(ioff2+kn)=a(ioff2+kn)+2.0_mpd*((2.0_mpd*
vecvk(l)*vtav-av(kn))*
vecvk(l)-av(kn)*
vecvk(l))
682 ioff2=roff(i+ioffp)+ioffp
684 a(ioff2+ifirst:ioff2+ilast)=a(ioff2+ifirst:ioff2+ilast)-2.0_mpd*
vecn(ifirst:ilast)*av(i)
685 a(ioff2+kn)=a(ioff2+kn)-2.0_mpd*
vecvk(l)*av(i)
693#ifdef SCOREP_USER_ENABLE
694 scorep_user_region_by_name_end(
"UR_qlssq")
719 INTEGER(mpi) :: ifirst
720 INTEGER(mpi) :: ilast
721 INTEGER(mpi) :: ifirst2
722 INTEGER(mpi) :: ilast2
723 INTEGER(mpl) :: ioff1
724 INTEGER(mpl) :: ioff2
725 INTEGER(mpi) :: istat(3)
735 INTEGER(mpl) :: length
741 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecAv
742 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matvtvp
743 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matvtAvp
744 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matCoeff
745 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: irangeCoeff
747 INTEGER(mpi),
INTENT(IN) :: m
748 REAL(mpd),
INTENT(IN OUT) :: B(npar*m-(m*m-m)/2)
749 LOGICAL,
INTENT(IN) :: t
752 SUBROUTINE aprod(n,l,x,is,ie,y)
754 INTEGER(mpi),
INTENT(in) :: n
755 INTEGER(mpl),
INTENT(in) :: l
756 REAL(mpd),
INTENT(IN) :: x(n)
757 INTEGER(mpi),
INTENT(in) :: is
758 INTEGER(mpi),
INTENT(in) :: ie
759 REAL(mpd),
INTENT(OUT) :: y(n)
763#ifdef SCOREP_USER_ENABLE
764 scorep_user_region_by_name_begin(
"UR_qlpssq", scorep_user_region_type_common)
768 CALL mpalloc(vecav,length,
'qlpssq: A*v')
770 CALL mpalloc(matvtvp,length,
"qlpssq: v^t*v'")
772 CALL mpalloc(matvtavp,length,
"qlpssq: v^t*(A*v')")
774 CALL mpalloc(matcoeff,length,
'qlpssq: coefficients')
777 CALL mpalloc(irangecoeff,2_mpl,length,
'qlpssq: non zero coefficient range')
787 ioff1=int(k-1,mpl)*int(
ncon,mpl)
788 irangecoeff(1,k)=
ncon
799 CALL aprod(npar,0_mpl,
vecn(1:npar),ifirst,ilast,vecav(1:npar))
800 CALL aprod(npar,0_mpl,
vecn(1:npar),kn,kn,vecav(1:npar))
806 ioff2=int(k2-1,mpl)*int(
ncon,mpl)
811 IF (kn >= ifirst2.AND.kn <= ilast2) v2kn=
matv(
ioffrow(k2)+1+kn-ifirst2)
813 l1=max(ifirst,ifirst2)
816 IF (l1 <= l2) vtvp=vtvp+dot_product(
vecn(l1:l2), &
819 IF (abs(vtvp) > 16.0_mpd*epsilon(vtvp))
THEN
820 matvtvp(ioff1+k2)=vtvp
821 matvtvp(ioff2+k)=vtvp
824 matvtvp(ioff1+k)=1.0_mpd
834 matvtavp(ioff1+k2)=
vecvk(k2)*vecav(kn2)+dot_product(vecav(ifirst2:ilast2), &
839 vtav=matvtavp(ioff1+k)
844 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*
vecn(i)*vtav-vecav(i))*
vecn(l)-vecav(l)*
vecn(i))
852 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*
vecn(l)
862 ioff1=int(k-1,mpl)*int(
ncon,mpl)
875 vtav=matvtavp(ioff1+k)+dot_product(matcoeff(ioff1+l1:ioff1+l2),matvtvp(ioff1+l1:ioff1+l2))
880 IF (matcoeff(ioff1+k2) == 0._mpd) cycle
881 if (istat(1)==0) istat(1)=k2
889 vecav(ifirst2:ilast2)=vecav(ifirst2:ilast2)+matcoeff(ioff1+k2)* &
891 vecav(kn2)=vecav(kn2)+matcoeff(ioff1+k2)*
vecvk(k2)
899 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*
vecn(i)*vtav-vecav(i))*
vecn(l)-vecav(l)*
vecn(i))
907 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*
vecn(l)
915 ioff2=int(k2-1,mpl)*int(
ncon,mpl)
916 vtvp=matvtvp(ioff1+k2)
920 vtavp=matvtavp(ioff2+k)
921 IF (l1 <= l2) vtavp=vtavp+dot_product(matcoeff(ioff2+l1:ioff2+l2),matvtvp(ioff1+l1:ioff1+l2))
924 matcoeff(ioff2+k)=matcoeff(ioff2+k)+2.0_mpd*(2.0_mpd*vtav*vtvp-vtavp)
925 IF (vtvp /= 0._mpd)
THEN
926 l1=min(l1,irangecoeff(1,k))
927 l2=max(l2,irangecoeff(2,k))
928 matcoeff(ioff2+l1:ioff2+l2)=matcoeff(ioff2+l1:ioff2+l2)-2.0_mpd*matcoeff(ioff1+l1:ioff1+l2)*vtvp
940#ifdef SCOREP_USER_ENABLE
941 scorep_user_region_by_name_end(
"UR_qlpssq")
960 INTEGER(mpi) :: ibpar
961 INTEGER(mpi) :: icoff
962 INTEGER(mpi) :: iclast
963 INTEGER(mpl) :: idiag
965 REAL(mpd),
INTENT(OUT) :: emin
966 REAL(mpd),
INTENT(OUT) :: emax
973 idiag=int(
ncon,mpl)*int(icoff,mpl)+1
975 IF (abs(emax) < abs(
matl(idiag))) emax=
matl(idiag)
976 IF (abs(emin) > abs(
matl(idiag))) emin=
matl(idiag)
995 INTEGER(mpi) :: icoff
996 INTEGER(mpi) :: iclast
997 INTEGER(mpl) :: idiag
999 INTEGER(mpi) :: nconb
1001 REAL(mpd),
INTENT(IN) :: d(ncon)
1002 REAL(mpd),
INTENT(OUT) :: y(ncon)
1008 idiag=int(ncon,mpl)*int(iclast-1,mpl)+nconb
1010 y(k)=(d(k)-dot_product(
matl(idiag+1:idiag+nconb-k),y(k+1:nconb)))/
matl(idiag)
1022 INTEGER(mpi),
INTENT(IN) :: ib
1035 INTEGER(mpi) :: ifirst
1036 INTEGER(mpi) :: ilast
1037 INTEGER(mpl) :: ioff1
1038 INTEGER(mpl) :: ioff2
1039 INTEGER(mpi) :: istat(6)
1054 v1=0.;v2=0.;v3=0.;v4=0.
1062 IF (
vecn(j) /= 0.0_mpd)
THEN
1064 IF (istat(3) == 0)
THEN
1074 IF (
matl(ioff2+j) /= 0.0_mpd)
THEN
1076 IF (istat(6) == 0)
THEN
1088100
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