SND@LHC Software
Loading...
Searching...
No Matches
mpqldec.f90
Go to the documentation of this file.
1
25
27MODULE mpqldec
28 USE mpdef
29 IMPLICIT NONE
30
31 INTEGER(mpi) :: npar
32 INTEGER(mpi) :: ncon
33 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matv
34 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matl
35 REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecn
36
37END MODULE mpqldec
38
44SUBROUTINE qlini(n,m)
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')
63END SUBROUTINE qlini
64
65! 141217 C. Kleinwort, DESY-FH1
83SUBROUTINE qldec(a)
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
140END SUBROUTINE qldec
141
142! 190312 C. Kleinwort, DESY-BELLE
162SUBROUTINE qldecb(a,nb,b)
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
271END SUBROUTINE qldecb
272
273
282SUBROUTINE qlmlq(x,m,t)
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
315END SUBROUTINE qlmlq
316
317
326SUBROUTINE qlmrq(x,m,t)
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
358END SUBROUTINE qlmrq
359
360
368SUBROUTINE qlsmq(x,t)
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
406END SUBROUTINE qlsmq
407
408
419SUBROUTINE qlssq(aprod,A,t)
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
483END SUBROUTINE qlssq
484
485
497SUBROUTINE qlpssq(aprod,B,m,t)
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
591END SUBROUTINE qlpssq
592
593
601SUBROUTINE qlgete(emin,emax)
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
620END SUBROUTINE qlgete
621
622
630SUBROUTINE qlbsub(d,y)
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
647END SUBROUTINE qlbsub
mpdalc::mpalloc
allocate array
Definition mpdalc.f90:36
mpdalc::mpdealloc
deallocate array
Definition mpdalc.f90:41
qlmrq
subroutine qlmrq(x, m, t)
Multiply right by Q(t).
Definition mpqldec.f90:327
qldecb
subroutine qldecb(a, nb, b)
QL decomposition (for disjoint block matrix).
Definition mpqldec.f90:163
qlsmq
subroutine qlsmq(x, t)
Similarity transformation by Q(t).
Definition mpqldec.f90:369
qldec
subroutine qldec(a)
QL decomposition.
Definition mpqldec.f90:84
qlmlq
subroutine qlmlq(x, m, t)
Multiply left by Q(t).
Definition mpqldec.f90:283
qlbsub
subroutine qlbsub(d, y)
Backward substitution.
Definition mpqldec.f90:631
qlgete
subroutine qlgete(emin, emax)
Get eigenvalues.
Definition mpqldec.f90:602
qlssq
subroutine qlssq(aprod, a, t)
Similarity transformation by Q(t).
Definition mpqldec.f90:420
qlpssq
subroutine qlpssq(aprod, b, m, t)
Partial similarity transformation by Q(t).
Definition mpqldec.f90:498
qlini
subroutine qlini(n, m)
Initialize QL decomposition.
Definition mpqldec.f90:45
mpdalc
(De)Allocate vectors and arrays.
Definition mpdalc.f90:24
mpdef
Definition of constants.
Definition mpdef.f90:24
mpdef::mpd
integer, parameter mpd
Definition mpdef.f90:38
mpqldec
QL data.
Definition mpqldec.f90:27
mpqldec::ncon
integer(mpi) ncon
number of constraints
Definition mpqldec.f90:32
mpqldec::matv
real(mpd), dimension(:), allocatable matv
unit normals (v_i) of Householder reflectors
Definition mpqldec.f90:33
mpqldec::matl
real(mpd), dimension(:), allocatable matl
lower diagonal matrix L
Definition mpqldec.f90:34
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