SND@LHC Software
Loading...
Searching...
No Matches
Modules | Functions/Subroutines | Variables
mpqldec.f90 File Reference

QL decompostion. More...

Go to the source code of this file.

Modules

module  mpqldec
 QL data.
 

Functions/Subroutines

subroutine qlini (n, m)
 Initialize QL decomposition.
 
subroutine qldec (a)
 QL decomposition.
 
subroutine qldecb (a, nb, b)
 QL decomposition (for disjoint block matrix).
 
subroutine qlmlq (x, m, t)
 Multiply left by Q(t).
 
subroutine qlmrq (x, m, t)
 Multiply right by Q(t).
 
subroutine qlsmq (x, t)
 Similarity transformation by Q(t).
 
subroutine qlssq (aprod, a, t)
 Similarity transformation by Q(t).
 
subroutine qlpssq (aprod, b, m, t)
 Partial similarity transformation by Q(t).
 
subroutine qlgete (emin, emax)
 Get eigenvalues.
 
subroutine qlbsub (d, y)
 Backward substitution.
 

Variables

integer(mpi) mpqldec::npar
 number of parameters
 
integer(mpi) mpqldec::ncon
 number of constraints
 
real(mpd), dimension(:), allocatable mpqldec::matv
 unit normals (v_i) of Householder reflectors
 
real(mpd), dimension(:), allocatable mpqldec::matl
 lower diagonal matrix L
 
real(mpd), dimension(:), allocatable mpqldec::vecn
 normal vector
 

Detailed Description

QL decompostion.

Author
Claus Kleinwort, DESY, 2015 (Claus.nosp@m..Kle.nosp@m.inwor.nosp@m.t@de.nosp@m.sy.de)
Copyright
Copyright (c) 2015 Deutsches Elektronen-Synchroton, Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY

This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details.

You should have received a copy of the GNU Library General Public License along with this program (see the file COPYING.LIB for more details); if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

QL decomposition of constraints matrix by Householder transformations for solution by elimination.

Definition in file mpqldec.f90.

Function/Subroutine Documentation

◆ qlbsub()

subroutine qlbsub ( real(mpd), dimension(ncon), intent(in)  d,
real(mpd), dimension(ncon), intent(out)  y 
)

Backward substitution.

Get y from L^t*y=d.

Parameters
[in]dNcon vector, resdiduals
[out]yNcon vector, solution

Definition at line 630 of file mpqldec.f90.

631 USE mpqldec
632
633 IMPLICIT NONE
634 INTEGER(mpl) :: idiag
635 INTEGER(mpi) :: k
636
637 REAL(mpd), INTENT(IN) :: d(ncon)
638 REAL(mpd), INTENT(OUT) :: y(ncon)
639
640 ! solve L*y=d by forward substitution
641 idiag=ncon*ncon
642 DO k=ncon,1,-1
643 y(k)=(d(k)-dot_product(matl(idiag+1:idiag+ncon-k),y(k+1:ncon)))/matl(idiag)
644 idiag=idiag-ncon-1
645 END DO
646
mpqldec
QL data.
Definition mpqldec.f90:27
mpqldec::ncon
integer(mpi) ncon
number of constraints
Definition mpqldec.f90:32
mpqldec::matl
real(mpd), dimension(:), allocatable matl
lower diagonal matrix L
Definition mpqldec.f90:34

◆ qldec()

subroutine qldec ( real(mpd), dimension(*), intent(in)  a)

QL decomposition.

QL decomposition with Householder transformations. Decompose N-By-M matrix A into orthogonal N-by-N matrix Q and a N-by-M matrix containing zeros except for a lower triangular M-by-M matrix L (at the bottom): 

         | 0 |
 A = Q * |   |
         | L |

The decomposition is stored in a N-by-M matrix matV containing the unit normal vectors v_i of the hyperplanes (Householder reflectors) defining Q. The lower triangular matrix L is stored in the M-by-M matrix matL.

Parameters
[in]aNpar-by-Ncon matrix

Definition at line 83 of file mpqldec.f90.

84 USE mpqldec
85 USE mpdalc
86
87 ! cost[dot ops] ~= Npar*Ncon*Ncon
88
89 IMPLICIT NONE
90 INTEGER(mpi) :: i
91 INTEGER(mpl) :: ioff1
92 INTEGER(mpl) :: ioff2
93 INTEGER(mpl) :: ioff3
94 INTEGER(mpi) :: k
95 INTEGER(mpi) :: kn
96 INTEGER(mpl) :: length
97 REAL(mpd) :: nrm
98 REAL(mpd) :: sp
99
100 REAL(mpd), INTENT(IN) :: a(*)
101
102 ! prepare
103 length=npar*ncon
104 matv=a(1:length)
105 matl=0.0_mpd
106
107 ! Householder procedure
108 DO k=ncon,1,-1
109 kn=npar+k-ncon
110 ! column offset
111 ioff1=(k-1)*npar
112 ! get column
113 vecn(1:kn)=matv(ioff1+1:ioff1+kn)
114 nrm = sqrt(dot_product(vecn(1:kn),vecn(1:kn)))
115 IF (nrm == 0.0_mpd) cycle
116 !
117 IF (vecn(kn) >= 0.0_mpd) THEN
118 vecn(kn)=vecn(kn)+nrm
119 ELSE
120 vecn(kn)=vecn(kn)-nrm
121 END IF
122 ! create normal vector
123 nrm = sqrt(dot_product(vecn(1:kn),vecn(1:kn)))
124 vecn(1:kn)=vecn(1:kn)/nrm
125 ! transformation
126 ioff2=0
127 DO i=1,k
128 sp=dot_product(vecn(1:kn),matv(ioff2+1:ioff2+kn))
129 matv(ioff2+1:ioff2+kn)=matv(ioff2+1:ioff2+kn)-2.0_mpd*vecn(1:kn)*sp
130 ioff2=ioff2+npar
131 END DO
132 ! store column of L
133 ioff3=(k-1)*ncon
134 matl(ioff3+k:ioff3+ncon)=matv(ioff1+kn:ioff1+npar)
135 ! store normal vector
136 matv(ioff1+1:ioff1+kn)=vecn(1:kn)
137 matv(ioff1+kn+1:ioff1+npar)=0.0_mpd
138 END DO
139
mpdalc
(De)Allocate vectors and arrays.
Definition mpdalc.f90:24
mpqldec::matv
real(mpd), dimension(:), allocatable matv
unit normals (v_i) of Householder reflectors
Definition mpqldec.f90:33
mpqldec::vecn
real(mpd), dimension(:), allocatable vecn
normal vector
Definition mpqldec.f90:35
mpqldec::npar
integer(mpi) npar
number of parameters
Definition mpqldec.f90:31

◆ qldecb()

subroutine qldecb ( real(mpd), dimension(*), intent(in)  a,
integer(mpi), intent(in)  nb,
integer(mpi), dimension(3,*), intent(in)  b 
)

QL decomposition (for disjoint block matrix).

QL decomposition with Householder transformations. Decompose N-By-M matrix A into orthogonal N-by-N matrix Q and a N-by-M matrix containing zeros except for a lower triangular M-by-M matrix L (at the bottom): 

         | 0 |
 A = Q * |   |
         | L |

The decomposition is stored in a N-by-M matrix matV containing the unit normal vectors v_i of the hyperplanes (Householder reflectors) defining Q. The lower triangular matrix L is stored in the M-by-M matrix matL.

Parameters
[in]ablock compressed Npar-by-Ncon matrix
[in]nbnumber of blocks
[in]b3-by-Ncon+1 matrix (with block definition)

Definition at line 162 of file mpqldec.f90.

163 USE mpqldec
164 USE mpdalc
165
166 ! cost[dot ops] ~= Npar*Ncon*Ncon
167
168 IMPLICIT NONE
169 INTEGER(mpi) :: i
170 INTEGER(mpi) :: ib
171 INTEGER(mpi) :: ifirst
172 INTEGER(mpi) :: ilast
173 INTEGER(mpl) :: ioff1
174 INTEGER(mpl) :: ioff2
175 INTEGER(mpl) :: ioff3
176 INTEGER(mpi) :: k
177 INTEGER(mpi) :: k1
178 INTEGER(mpi) :: kn
179 INTEGER(mpi) :: ncb
180 INTEGER(mpi) :: npb
181 INTEGER(mpl) :: length
182 REAL(mpd) :: nrm
183 REAL(mpd) :: sp
184
185 REAL(mpd), INTENT(IN) :: a(*)
186 INTEGER(mpi), INTENT(IN) :: nb
187 INTEGER(mpi), INTENT(IN) :: b(3,*)
188
189 ! prepare
190 length=npar*ncon
191 matv=0.0_mpd
192 matl=0.0_mpd
193 ! expand a into matV
194 ioff1=0
195 ioff2=0
196 DO ib=1,nb
197 ncb=b(1,ib+1)-b(1,ib) ! number of constraints in block
198 npb=b(3,ib)+1-b(2,ib) ! number of parameters in block
199 ifirst=b(2,ib)
200 ilast=b(3,ib)
201 DO i=1,ncb
202 matv(ioff1+ifirst:ioff1+ilast)=a(ioff2+1:ioff2+npb)
203 ioff1=ioff1+npar
204 ioff2=ioff2+npb
205 END DO
206 END DO
207
208 ib=nb ! start with last block
209 k1=b(1,ib) ! first constraint in block
210 ! Householder procedure
211 DO k=ncon,1,-1
212 kn=npar+k-ncon
213 ! different block?
214 IF (k < k1) THEN
215 ib=ib-1
216 k1=b(1,ib)
217 END IF
218 ! index if first non-zero element
219 ifirst=b(2,ib)
220 IF (ifirst > kn) cycle
221 ! index of last element
222 ilast=min(b(3,ib),kn)
223 ! column offsets
224 ioff1=(k-1)*npar
225 ioff2=(k1-1)*npar
226 ! get column
227 vecn(kn)=0.0_mpd
228 vecn(ifirst:ilast)=matv(ioff1+ifirst:ioff1+ilast)
229 nrm = sqrt(dot_product(vecn(ifirst:ilast),vecn(ifirst:ilast)))
230 IF (nrm == 0.0_mpd) cycle
231 !
232 IF (vecn(kn) >= 0.0_mpd) THEN
233 vecn(kn)=vecn(kn)+nrm
234 ELSE
235 vecn(kn)=vecn(kn)-nrm
236 END IF
237
238 IF (ilast < kn) THEN
239 ! create normal vector
240 nrm = sqrt(dot_product(vecn(ifirst:ilast),vecn(ifirst:ilast))+vecn(kn)*vecn(kn))
241 vecn(ifirst:ilast)=vecn(ifirst:ilast)/nrm
242 vecn(kn)=vecn(kn)/nrm
243 ! transformation
244 DO i=k1,k
245 sp=dot_product(vecn(ifirst:ilast),matv(ioff2+ifirst:ioff2+ilast))
246 matv(ioff2+ifirst:ioff2+ilast)=matv(ioff2+ifirst:ioff2+ilast)-2.0_mpd*vecn(ifirst:ilast)*sp
247 matv(ioff2+kn)=-2.0_mpd*vecn(kn)*sp
248 ioff2=ioff2+npar
249 END DO
250 ELSE
251 ! create normal vector
252 nrm = sqrt(dot_product(vecn(ifirst:ilast),vecn(ifirst:ilast)))
253 vecn(ifirst:ilast)=vecn(ifirst:ilast)/nrm
254 ! transformation
255 DO i=k1,k
256 sp=dot_product(vecn(ifirst:ilast),matv(ioff2+ifirst:ioff2+ilast))
257 matv(ioff2+ifirst:ioff2+ilast)=matv(ioff2+ifirst:ioff2+ilast)-2.0_mpd*vecn(ifirst:ilast)*sp
258 ioff2=ioff2+npar
259 END DO
260 END IF
261
262 ! store column of L
263 ioff3=(k-1)*ncon
264 matl(ioff3+k:ioff3+ncon)=matv(ioff1+kn:ioff1+npar)
265 ! store normal vector
266 matv(ioff1+1:ioff1+npar)=0.0_mpd
267 matv(ioff1+ifirst:ioff1+ilast)=vecn(ifirst:ilast)
268 matv(ioff1+kn)=vecn(kn)
269 END DO
270

◆ qlgete()

subroutine qlgete ( real(mpd), intent(out)  emin,
real(mpd), intent(out)  emax 
)

Get eigenvalues.

Get smallest and largest |eigenvalue| of L.

Parameters
[out]emineigenvalue with smallest absolute value
[out]emaxeigenvalue with largest absolute value

Definition at line 601 of file mpqldec.f90.

602 USE mpqldec
603
604 IMPLICIT NONE
605 INTEGER(mpi) :: i
606 INTEGER(mpl) :: idiag
607
608 REAL(mpd), INTENT(OUT) :: emin
609 REAL(mpd), INTENT(OUT) :: emax
610
611 idiag=1
612 emax=matl(1)
613 emin=emax
614 DO i=2,ncon
615 idiag=idiag+ncon+1
616 IF (abs(emax) < abs(matl(idiag))) emax=matl(idiag)
617 IF (abs(emin) > abs(matl(idiag))) emin=matl(idiag)
618 END DO
619

◆ qlini()

subroutine qlini ( integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m 
)

Initialize QL decomposition.

Parameters
[in]nnumber of rows (parameters)
[in]mnumber of columns (constraints)

Definition at line 44 of file mpqldec.f90.

45 USE mpqldec
46 USE mpdalc
47
48 IMPLICIT NONE
49 INTEGER(mpl) :: length
50
51 INTEGER(mpi), INTENT(IN) :: n
52 INTEGER(mpi), INTENT(IN) :: m
53
54 npar=n
55 ncon=m
56 ! allocate
57 length=npar*ncon
58 CALL mpalloc(matv,length,'QLDEC: V')
59 length=ncon*ncon
60 CALL mpalloc(matl,length,'QLDEC: L')
61 length=npar
62 CALL mpalloc(vecn,length,'QLDEC: v')
mpdalc::mpalloc
allocate array
Definition mpdalc.f90:36

◆ qlmlq()

subroutine qlmlq ( real(mpd), dimension(*), intent(inout)  x,
integer(mpi), intent(in)  m,
logical, intent(in)  t 
)

Multiply left by Q(t).

Multiply left by Q(t) from QL decomposition.

Parameters
[in,out]xNpar-by-M matrix, overwritten with Q*X (t=false) or Q^t*X (t=true)
[in]mnumber of columns
[in]tuse transposed of Q

Definition at line 282 of file mpqldec.f90.

283 USE mpqldec
284
285 ! cost[dot ops] ~= N*M*Nhr
286
287 IMPLICIT NONE
288 INTEGER(mpi) :: i
289 INTEGER(mpl) :: ioff1
290 INTEGER(mpl) :: ioff2
291 INTEGER(mpi) :: j
292 INTEGER(mpi) :: k
293 INTEGER(mpi) :: kn
294 REAL(mpd) :: sp
295
296 REAL(mpd), INTENT(IN OUT) :: x(*)
297 INTEGER(mpi), INTENT(IN) :: m
298 LOGICAL, INTENT(IN) :: t
299
300 DO j=1,ncon
301 k=j
302 IF (t) k=ncon+1-j
303 kn=npar+k-ncon
304 ! column offset
305 ioff1=(k-1)*npar
306 ! transformation
307 ioff2=0
308 DO i=1,m
309 sp=dot_product(matv(ioff1+1:ioff1+kn),x(ioff2+1:ioff2+kn))
310 x(ioff2+1:ioff2+kn)=x(ioff2+1:ioff2+kn)-2.0_mpd*matv(ioff1+1:ioff1+kn)*sp
311 ioff2=ioff2+npar
312 END DO
313 END DO
314

◆ qlmrq()

subroutine qlmrq ( real(mpd), dimension(*), intent(inout)  x,
integer(mpi), intent(in)  m,
logical, intent(in)  t 
)

Multiply right by Q(t).

Multiply right by Q(t) from QL decomposition.

Parameters
[in,out]xM-by-Npar matrix, overwritten with X*Q (t=false) or X*Q^t (t=true)
[in]mnumber of rows
[in]tuse transposed of Q

Definition at line 326 of file mpqldec.f90.

327 USE mpqldec
328
329 ! cost[dot ops] ~= N*M*Nhr
330
331 IMPLICIT NONE
332 INTEGER(mpi) :: i
333 INTEGER(mpl) :: ioff1
334 INTEGER(mpl) :: iend
335 INTEGER(mpi) :: j
336 INTEGER(mpi) :: k
337 INTEGER(mpi) :: kn
338 REAL(mpd) :: sp
339
340 REAL(mpd), INTENT(IN OUT) :: x(*)
341 INTEGER(mpi), INTENT(IN) :: m
342 LOGICAL, INTENT(IN) :: t
343
344 DO j=1,ncon
345 k=j
346 IF (.not.t) k=ncon+1-j
347 kn=npar+k-ncon
348 ! column offset
349 ioff1=(k-1)*npar
350 ! transformation
351 iend=m*kn
352 DO i=1,npar
353 sp=dot_product(matv(ioff1+1:ioff1+kn),x(i:iend:m))
354 x(i:iend:m)=x(i:iend:m)-2.0_mpd*matv(ioff1+1:ioff1+kn)*sp
355 END DO
356 END DO
357

◆ qlpssq()

subroutine qlpssq ( external subroutine(integer(mpi), intent(in) n, real(mpd), dimension(n), intent(in) x, real(mpd), dimension(n), intent(out) y)  aprod,
real(mpd), dimension(*), intent(inout)  b,
integer(mpi), intent(in)  m,
logical, intent(in)  t 
)

Partial similarity transformation by Q(t).

Partial similarity transformation for symmetric matrix by Q from QL decomposition. Calculate corrections to band part of matrix.

Parameters
[in]aprodexternal procedure to calculate A*v
[in,out]Bband part of symmetric Npar-by-Npar matrix A in symmetric storage mode, overwritten with band part of Q^t*A*Q (t=false) or Q^t*A*Q (t=true)
[in]mband width (including diagonal)
[in]tuse transposed of Q

Definition at line 497 of file mpqldec.f90.

498 USE mpqldec
499 USE mpdalc
500
501 ! cost[dot ops] ~= N*N*Nhr
502
503 IMPLICIT NONE
504 INTEGER(mpi) :: i
505 INTEGER(mpl) :: ioff1
506 INTEGER(mpl) :: ioff2
507 INTEGER(mpi) :: j
508 INTEGER(mpi) :: j2
509 INTEGER(mpi) :: k
510 INTEGER(mpi) :: k2
511 INTEGER(mpi) :: kn
512 INTEGER(mpi) :: l
513 INTEGER(mpl) :: length
514 INTEGER(mpi) :: mbnd
515 REAL(mpd) :: vtAv
516 REAL(mpd) :: vtAvp
517 REAL(mpd) :: vtvp
518 REAL(mpd), DIMENSION(:), ALLOCATABLE :: Av ! A*v
519
520 REAL(mpd), INTENT(IN OUT) :: B(*)
521 INTEGER(mpi), INTENT(IN) :: m
522 LOGICAL, INTENT(IN) :: t
523
524 INTERFACE
525 SUBROUTINE aprod(n,x,y) ! y=A*x
526 USE mpdef
527 INTEGER(mpi), INTENT(in) :: n
528 REAL(mpd), INTENT(IN) :: x(n)
529 REAL(mpd), INTENT(OUT) :: y(n)
530 END SUBROUTINE aprod
531 END INTERFACE
532
533 length=npar
534 length=npar*ncon
535 CALL mpalloc(av,length,'qlpssq: Av')
536
537 mbnd=max(0,m-1) ! band width without diagonal
538 ! A*V
539 ioff1=0
540 DO i=1,ncon
541 CALL aprod(npar,matv(ioff1+1:ioff1+npar),av(ioff1+1:ioff1+npar))
542 ioff1=ioff1+npar
543 END DO
544
545 DO j=1,ncon
546 k=j
547 IF (t) k=ncon+1-j
548 kn=npar+k-ncon
549 ! column offset
550 ioff1=(k-1)*npar
551 ! transformation (diagonal block)
552 ! diagonal block
553 ! v^t*A*v
554 vtav=dot_product(matv(ioff1+1:ioff1+kn),av(ioff1+1:ioff1+kn))
555 ! update
556 ioff2=0
557 DO i=1,kn
558 ! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
559 DO l=max(1,i-mbnd),i
560 ioff2=ioff2+1
561 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*matv(ioff1+i)*vtav-av(ioff1+i))*matv(ioff1+l)-av(ioff1+l)*matv(ioff1+i))
562 END DO
563 END DO
564 ! off diagonal block
565 DO i=kn+1,npar
566 ! correct with -2Av*v^t
567 DO l=max(1,i-mbnd),i
568 ioff2=ioff2+1
569 b(ioff2)=b(ioff2)-2.0_mpd*av(ioff1+i)*matv(ioff1+l)
570 END DO
571 END DO
572 ! correct A*v for the remainung v
573 DO j2=j+1,ncon
574 k2=j2
575 IF (t) k2=ncon+1-j2
576 ioff2=(k2-1)*npar
577 vtvp=dot_product(matv(ioff1+1:ioff1+npar),matv(ioff2+1:ioff2+npar)) ! v^t*v'
578 vtavp=dot_product(matv(ioff1+1:ioff1+npar),av(ioff2+1:ioff2+npar)) ! v^t*(A*v')
579 DO i=1,kn
580 av(ioff2+i)=av(ioff2+i)+2.0_mpd*((2.0_mpd*matv(ioff1+i)*vtav-av(ioff1+i))*vtvp-matv(ioff1+i)*vtavp)
581 END DO
582 DO i=kn+1,npar
583 av(ioff2+i)=av(ioff2+i)-2.0_mpd*av(ioff1+i)*vtvp
584 END DO
585 END DO
586
587 END DO
588
589 CALL mpdealloc(av)
590
mpdalc::mpdealloc
deallocate array
Definition mpdalc.f90:41
mpdef
Definition of constants.
Definition mpdef.f90:24

◆ qlsmq()

subroutine qlsmq ( real(mpd), dimension(*), intent(inout)  x,
logical, intent(in)  t 
)

Similarity transformation by Q(t).

Similarity transformation by Q from QL decomposition.

Parameters
[in,out]xNpar-by-Npar matrix, overwritten with Q*X*Q^t (t=false) or Q^t*X*Q (t=true)
[in]tuse transposed of Q

Definition at line 368 of file mpqldec.f90.

369 USE mpqldec
370
371 ! cost[dot ops] ~= N*N*Nhr
372
373 IMPLICIT NONE
374 INTEGER(mpi) :: i
375 INTEGER(mpl) :: ioff1
376 INTEGER(mpl) :: ioff2
377 INTEGER(mpl) :: iend
378 INTEGER(mpi) :: j
379 INTEGER(mpi) :: k
380 INTEGER(mpi) :: kn
381 REAL(mpd) :: sp
382
383 REAL(mpd), INTENT(IN OUT) :: x(*)
384 LOGICAL, INTENT(IN) :: t
385
386 DO j=1,ncon
387 k=j
388 IF (t) k=ncon+1-j
389 kn=npar+k-ncon
390 ! column offset
391 ioff1=(k-1)*npar
392 ! transformation
393 iend=npar*kn
394 DO i=1,npar
395 sp=dot_product(matv(ioff1+1:ioff1+kn),x(i:iend:npar))
396 x(i:iend:npar)=x(i:iend:npar)-2.0_mpd*matv(ioff1+1:ioff1+kn)*sp
397 END DO
398 ioff2=0
399 DO i=1,npar
400 sp=dot_product(matv(ioff1+1:ioff1+kn),x(ioff2+1:ioff2+kn))
401 x(ioff2+1:ioff2+kn)=x(ioff2+1:ioff2+kn)-2.0_mpd*matv(ioff1+1:ioff1+kn)*sp
402 ioff2=ioff2+npar
403 END DO
404 END DO
405

◆ qlssq()

subroutine qlssq ( external subroutine(integer(mpi), intent(in) n, real(mpd), dimension(n), intent(in) x, real(mpd), dimension(n), intent(out) y)  aprod,
real(mpd), dimension(*), intent(inout)  a,
logical, intent(in)  t 
)

Similarity transformation by Q(t).

Similarity transformation for symmetric matrix by Q from QL decomposition.

Parameters
[in]aprodexternal procedure to calculate A*v
[in,out]Asymmetric Npar-by-Npar matrix A in symmetric storage mode (V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, ...), overwritten with Q*A*Q^t (t=false) or Q^t*A*Q (t=true)
[in]tuse transposed of Q

Definition at line 419 of file mpqldec.f90.

420 USE mpqldec
421 USE mpdalc
422
423 ! cost[dot ops] ~= N*N*Nhr
424
425 IMPLICIT NONE
426 INTEGER(mpi) :: i
427 INTEGER(mpl) :: ioff1
428 INTEGER(mpl) :: ioff2
429 INTEGER(mpi) :: j
430 INTEGER(mpi) :: k
431 INTEGER(mpi) :: kn
432 INTEGER(mpi) :: l
433 INTEGER(mpl) :: length
434 REAL(mpd) :: vtAv
435 REAL(mpd), DIMENSION(:), ALLOCATABLE :: Av
436
437 REAL(mpd), INTENT(IN OUT) :: A(*)
438 LOGICAL, INTENT(IN) :: t
439
440 INTERFACE
441 SUBROUTINE aprod(n,x,y) ! y=A*x
442 USE mpdef
443 INTEGER(mpi), INTENT(in) :: n
444 REAL(mpd), INTENT(IN) :: x(n)
445 REAL(mpd), INTENT(OUT) :: y(n)
446 END SUBROUTINE aprod
447 END INTERFACE
448
449 length=npar
450 CALL mpalloc(av,length,'qlssq: A*v')
451
452 DO j=1,ncon
453 k=j
454 IF (t) k=ncon+1-j
455 kn=npar+k-ncon
456 ! column offset
457 ioff1=(k-1)*npar
458 ! A*v
459 CALL aprod(npar,matv(ioff1+1:ioff1+npar),av(1:npar))
460 ! transformation
461 ! diagonal block
462 ! v^t*A*v
463 vtav=dot_product(matv(ioff1+1:ioff1+kn),av(1:kn))
464 ! update
465 ioff2=0
466 DO i=1,kn
467 ! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
468 DO l=1,i
469 ioff2=ioff2+1
470 a(ioff2)=a(ioff2)+2.0_mpd*((2.0_mpd*matv(ioff1+i)*vtav-av(i))*matv(ioff1+l)-av(l)*matv(ioff1+i))
471 END DO
472 END DO
473 ! off diagonal block
474 DO i=kn+1,npar
475 ! correct with -2Av*v^t
476 a(ioff2+1:ioff2+kn)=a(ioff2+1:ioff2+kn)-2.0_mpd*matv(ioff1+1:ioff1+kn)*av(i)
477 ioff2=ioff2+i
478 END DO
479 END DO
480
481 CALL mpdealloc(av)
482