Millepede-II V04-16-00
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 INTEGER(mpi) :: nblock
34 INTEGER(mpl) :: matsize
35 INTEGER(mpi) :: iblock
36 INTEGER(mpi) :: monpg
37 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matv
38 REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecvk
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
45 INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: irangeparnz
46
47END MODULE mpqldec
48
57SUBROUTINE qlini(n,m,l,s,k)
58 USE mpqldec
59 USE mpdalc
60
61 IMPLICIT NONE
62 INTEGER(mpl) :: length
63
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
69
70 npar=n
71 ncon=m
72 nblock=l
73 matsize=s
74 iblock=1
75 monpg=k
76 ! allocate
77 length=matsize
78 !print *, ' full length (V)', length
79 CALL mpalloc(matv,length,'QLDEC: V')
80 matv=0.
81 length=int(ncon,mpl)*int(ncon,mpl)
82 CALL mpalloc(matl,length,'QLDEC: L')
83 matl=0.
84 length=npar
85 CALL mpalloc(vecn,length,'QLDEC: v')
86 length=ncon
87 CALL mpalloc(vecvk,length,'QLDEC: sec. diag(V)')
88 vecvk=0.
89 CALL mpalloc(ioffpar,length,'QLDEC: parameter offsets (V)')
90 ioffpar=0
91 CALL mpalloc(irangeparnz,2_mpl,length,'QLDEC: parameter non zero range (V)')
92 length=ncon+1
93 CALL mpalloc(ioffrow,length,'QLDEC: row offsets (V)')
94 ioffrow=0
95 length=nblock
96 CALL mpalloc(nparblock,length,'QLDEC: npar in block')
97 nparblock=0
98 length=nblock+1
99 CALL mpalloc(ioffblock,length,'QLDEC: ioff for block')
100 ioffblock=0
101END SUBROUTINE qlini
102
103! 141217 C. Kleinwort, DESY-FH1
121SUBROUTINE qldec(a)
122 USE mpqldec
123 USE mpdalc
124
125 ! cost[dot ops] ~= Npar*Ncon*Ncon
126
127 IMPLICIT NONE
128 INTEGER(mpi) :: i
129 INTEGER(mpl) :: ioff1
130 INTEGER(mpl) :: ioff2
131 INTEGER(mpl) :: ioff3
132 INTEGER(mpi) :: k
133 INTEGER(mpi) :: kn
134 INTEGER(mpl) :: length
135 REAL(mpd) :: nrm
136 REAL(mpd) :: sp
137
138 REAL(mpd), INTENT(IN) :: a(matsize)
139
140 ! prepare
141 vecvk=0._mpd
142 length=int(npar,mpl)*int(ncon,mpl)
143 matv=a(1:length)
144 matl=0.0_mpd
145 ! implemented as single block
146 nblock=1
147 nparblock(1)=npar
148 ioffblock(2)=ncon
149 DO k=1,ncon
150 ioffrow(k+1)=ioffrow(k)+npar
151 END DO
152
153 ! Householder procedure
154 DO k=ncon,1,-1
155 ! monitoring
156 IF(monpg>0) CALL monpgs(ncon+1-k)
157 kn=npar+k-ncon
158 ! column offset
159 ioff1=int(k-1,mpl)*int(npar,mpl)
160 ! get column
161 vecn(1:kn)=matv(ioff1+1:ioff1+kn)
162 nrm = sqrt(dot_product(vecn(1:kn),vecn(1:kn)))
163 IF (nrm == 0.0_mpd) cycle
164 !
165 IF (vecn(kn) >= 0.0_mpd) THEN
166 vecn(kn)=vecn(kn)+nrm
167 ELSE
168 vecn(kn)=vecn(kn)-nrm
169 END IF
170 ! create normal vector
171 nrm = sqrt(dot_product(vecn(1:kn),vecn(1:kn)))
172 vecn(1:kn)=vecn(1:kn)/nrm
173 ! transformation
174 ioff2=0
175 DO i=1,k
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
178 ioff2=ioff2+npar
179 END DO
180 ! store column of L
181 ioff3=int(k-1,mpl)*int(ncon,mpl)
182 matl(ioff3+k:ioff3+ncon)=matv(ioff1+kn:ioff1+npar)
183 ! store normal vector
184 matv(ioff1+1:ioff1+kn-1)=vecn(1:kn-1)
185 matv(ioff1+kn:ioff1+npar)=0.0_mpd
186 irangeparnz(1,k)=1
187 irangeparnz(2,k)=kn-1
188 ! store secondary diagonal
189 vecvk(k)=vecn(kn)
190 END DO
191
192END SUBROUTINE qldec
193
194! 190312 C. Kleinwort, DESY-BELLE
215SUBROUTINE qldecb(a,bpar,bcon,rcon)
216 USE mpqldec
217 USE mpdalc
218
219 ! cost[dot ops] ~= Npar*Ncon*Ncon
220
221 IMPLICIT NONE
222 INTEGER(mpi) :: i
223 INTEGER(mpi) :: ibcon
224 INTEGER(mpi) :: iblast
225 INTEGER(mpi) :: iboff
226 INTEGER(mpi) :: ibpar
227 INTEGER(mpi) :: ifirst
228 INTEGER(mpi) :: ilast
229 INTEGER(mpi) :: in
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
237 INTEGER(mpi) :: k
238 INTEGER(mpi) :: k1
239 INTEGER(mpi) :: kn
240 INTEGER(mpi) :: ncb
241 INTEGER(mpi) :: ncol
242 INTEGER(mpi) :: npb
243 REAL(mpd) :: nrm
244 REAL(mpd) :: sp
245
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)
250
251 !$POMP INST BEGIN(qldecb)
252 ! prepare
253 vecvk=0.0_mpd
254 matv=a(1:matsize)
255 matl=0.0_mpd
256
257 ! prepare offsets
258 icoff=0
259 DO ibpar=1,nblock ! parameter block
260 iclast=icoff
261 DO ibcon=bpar(2,ibpar)+1, bpar(2,ibpar+1)! constraint block
262 ncb=bcon(1,ibcon+1)-bcon(1,ibcon) ! number of constraints in constraint block
263 npb=bcon(3,ibcon)+1-bcon(2,ibcon) ! number of parameters in constraint block
264 ifirst=bcon(2,ibcon)
265 ilast=bcon(3,ibcon)
266 DO i=bcon(1,ibcon),bcon(1,ibcon+1)-1
267 ! non-zero range: first, last parameter
268 irangeparnz(1,i)=rcon(1,i)
269 irangeparnz(2,i)=rcon(2,i)
270 ! storage: parameter, row offset
271 ioffpar(i)=rcon(3,i)-1
272 ioffrow(i+1)=ioffrow(i)+rcon(4,i)-ioffpar(i)
273 END DO
274 iclast=iclast+ncb
275 END DO
276 ! set up matL
277 iblast=bpar(1,ibpar+1) ! last parameter in parameter block
278 DO k=icoff+1,iclast
279 kn=iblast+k-iclast
280 ioff1=ioffrow(k)
281 npb=int(ioffrow(k+1)-ioff1,mpi)
282 ! size of overlap
283 ncol=ioffpar(k)+npb-kn
284 IF (ncol >= 0) THEN
285 ioff3=int(k-1,mpl)*int(ncon,mpl)
286 matl(ioff3+iclast-ncol-icoff:ioff3+iclast-icoff)=matv(ioff1+npb-ncol:ioff1+npb)
287 END IF
288 END DO
289 icoff=iclast
290 nparblock(ibpar)=bpar(1,ibpar+1)-bpar(1,ibpar)
291 ioffblock(ibpar+1)=icoff
292 END DO
293
294 DO ibpar=1,nblock ! parameter block
295 iboff=bpar(1,ibpar) ! last offset in parameter block
296 iblast=bpar(1,ibpar+1) ! last parameter in parameter block
297 icoff=ioffblock(ibpar) ! constraint offset in parameter block
298 iclast=ioffblock(ibpar+1) ! last constraint in parameter block
299 IF(iclast <= icoff) cycle ! no constraints
300 ibcon=bpar(2,ibpar+1) ! start with last constraint block
301 k1=bcon(1,ibcon) ! first constraint in block
302 ! Householder procedure
303 DO k=iclast,icoff+1,-1
304 ! monitoring
305 IF(monpg>0) CALL monpgs(ncon+1-k)
306 kn=iblast+k-iclast
307 ! different constraint block?
308 IF (k < k1) THEN
309 ibcon=ibcon-1
310 k1=bcon(1,ibcon)
311 END IF
312 ! parameter offset
313 ipoff=ioffpar(k)
314 ! index if first non-zero parameter
315 ifirst=ipoff+1
316 IF (ifirst > kn) cycle
317 ! column offset
318 ioff1=ioffrow(k)
319 ! number of parameter
320 npb=int(ioffrow(k+1)-ioff1,mpi)
321 ! index of last parameter
322 iplast=ioffpar(k)+npb
323 ! index of last used parameter
324 ilast=min(iplast,kn)
325 ! number of used columns
326 ncol=ilast-ipoff
327 ! get column
328 in=iblast+k1-iclast
329 vecn(in:kn)=0.0_mpd
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
333 !
334 IF (vecn(kn) >= 0.0_mpd) THEN
335 vecn(kn)=vecn(kn)+nrm
336 ELSE
337 vecn(kn)=vecn(kn)-nrm
338 END IF
339 ! create normal vector
340 IF (ilast < kn) THEN
341 nrm = sqrt(dot_product(vecn(ifirst:ilast),vecn(ifirst:ilast))+vecn(kn)*vecn(kn))
342 vecn(ifirst:ilast)=vecn(ifirst:ilast)/nrm
343 vecn(kn)=vecn(kn)/nrm
344 ELSE
345 nrm = sqrt(dot_product(vecn(ifirst:ilast),vecn(ifirst:ilast)))
346 vecn(ifirst:ilast)=vecn(ifirst:ilast)/nrm
347 END IF
348 ! update L too
349 ioff3=int(k1-2,mpl)*int(ncon,mpl)
350 ! transformation
351 DO i=k1,k
352 ioff3=ioff3+ncon
353 IF (irangeparnz(1,k) > irangeparnz(2,i)) cycle ! no overlap
354 ioff2=ioffrow(i)+ioffpar(k)-ioffpar(i)
355 sp=dot_product(vecn(ifirst:ilast),matv(ioff2+1:ioff2+ncol))
356 IF (sp /= 0.0_mpd) THEN
357 ! update matV
358 matv(ioff2+1:ioff2+ncol)=matv(ioff2+1:ioff2+ncol)-2.0_mpd*vecn(ifirst:ilast)*sp
359 ! update matL
360 in=iblast+i-iclast
361 matl(ioff3+i-icoff:ioff3+k-icoff)=matl(ioff3+i-icoff:ioff3+k-icoff)-2.0_mpd*vecn(in:kn)*sp
362 ! update non zero range
363 irangeparnz(1,i)=min(irangeparnz(1,i),irangeparnz(1,k))
364 irangeparnz(2,i)=max(irangeparnz(2,i),irangeparnz(2,k))
365 END IF
366 END DO
367 ! store secondary diagonal
368 vecvk(icoff+k)=vecn(kn)
369 ! store normal vector (non zero part)
370 ifirst=irangeparnz(1,k)
371 ilast=min(irangeparnz(2,k),kn-1)
372 ncol=ilast-ifirst+1
373 matv(ioff1+1:ioff1+ncol)=vecn(ifirst:ilast)
374 matv(ioff1+ncol+1:ioff1+npb)=0.0_mpd
375 ! local to parameter block
376 irangeparnz(1,k)=ifirst-iboff
377 irangeparnz(2,k)=ilast-iboff
378 END DO
379 END DO
380 !$POMP INST END(qldecb)
381
382
383END SUBROUTINE qldecb
384
385
394SUBROUTINE qlmlq(x,m,t)
395 USE mpqldec
396
397 ! cost[dot ops] ~= N*M*Nhr
398
399 IMPLICIT NONE
400 INTEGER(mpi) :: i
401 INTEGER(mpi) :: icoff
402 INTEGER(mpi) :: iclast
403 INTEGER(mpi) :: ifirst
404 INTEGER(mpi) :: ilast
405 INTEGER(mpl) :: ioff2
406 INTEGER(mpi) :: j
407 INTEGER(mpi) :: k
408 INTEGER(mpi) :: l
409 INTEGER(mpi) :: kn
410 INTEGER(mpi) :: nconb
411 INTEGER(mpi) :: nparb
412 REAL(mpd) :: sp
413
414 INTEGER(mpi), INTENT(IN) :: m
415 REAL(mpd), INTENT(IN OUT) :: x(INT(npar,mpl)*INT(m,mpl))
416 LOGICAL, INTENT(IN) :: t
417
418 icoff=ioffblock(iblock) ! constraint offset in parameter block
419 iclast=ioffblock(iblock+1) ! last constraint in parameter block
420 nconb=iclast-icoff ! number of constraints in block
421 nparb=nparblock(iblock) ! number of parameters in block
422 DO j=1,nconb
423 k=j
424 IF (t) k=nconb+1-j
425 kn=nparb+k-nconb
426 ! expand row 'l' of matV into vecN
427 l=k+icoff
428 ! non-zero range (excluding 'kn')
429 ifirst=irangeparnz(1,l)
430 ilast=irangeparnz(2,l)
431 vecn(1:nparb)=0._mpd
432 vecn(ifirst:ilast)=matv(ioffrow(l)+1:ioffrow(l)+1+ilast-ifirst)
433 vecn(kn)=vecvk(l)
434 ! transformation
435 ioff2=0
436 DO i=1,m
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
440 ioff2=ioff2+nparb
441 END DO
442 END DO
443
444END SUBROUTINE qlmlq
445
446
455SUBROUTINE qlmrq(x,m,t)
456 USE mpqldec
457
458 ! cost[dot ops] ~= N*M*Nhr
459
460 IMPLICIT NONE
461 INTEGER(mpi) :: i
462 INTEGER(mpl) :: iend
463 INTEGER(mpi) :: ifirst
464 INTEGER(mpi) :: ilast
465 INTEGER(mpi) :: j
466 INTEGER(mpi) :: k
467 INTEGER(mpi) :: kn
468 REAL(mpd) :: sp
469
470 INTEGER(mpi), INTENT(IN) :: m
471 REAL(mpd), INTENT(IN OUT) :: x(INT(m,mpl)*INT(npar,mpl))
472 LOGICAL, INTENT(IN) :: t
473
474 DO j=1,ncon
475 k=j
476 IF (.not.t) k=ncon+1-j
477 kn=npar+k-ncon
478 ! expand row 'k' of matV into vecN
479 ! non-zero range (excluding 'kn')
480 ifirst=irangeparnz(1,k)
481 ilast=irangeparnz(2,k)
482 vecn=0._mpd
483 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
484 vecn(kn)=vecvk(k)
485 ! transformation
486 iend=m*kn
487 DO i=1,npar
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
490 END DO
491 END DO
492
493END SUBROUTINE qlmrq
494
495
503SUBROUTINE qlsmq(x,t)
504 USE mpqldec
505
506 ! cost[dot ops] ~= N*N*Nhr
507
508 IMPLICIT NONE
509 INTEGER(mpi) :: i
510 INTEGER(mpl) :: ioff2
511 INTEGER(mpl) :: iend
512 INTEGER(mpi) :: ifirst
513 INTEGER(mpi) :: ilast
514 INTEGER(mpi) :: j
515 INTEGER(mpi) :: k
516 INTEGER(mpi) :: kn
517 REAL(mpd) :: sp
518
519 REAL(mpd), INTENT(IN OUT) :: x(INT(npar,mpl)*INT(npar,mpl))
520 LOGICAL, INTENT(IN) :: t
521
522 DO j=1,ncon
523 ! monitoring
524 IF(monpg>0) CALL monpgs(j)
525 k=j
526 IF (t) k=ncon+1-j
527 kn=npar+k-ncon
528 ! expand row 'k' of matV into vecN
529 ! non-zero range (excluding 'kn')
530 ifirst=irangeparnz(1,k)
531 ilast=irangeparnz(2,k)
532 vecn=0._mpd
533 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
534 vecn(kn)=vecvk(k)
535 ! transformation
536 iend=int(npar,mpl)*int(kn,mpl)
537 DO i=1,npar
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
540 END DO
541 ioff2=0
542 DO i=1,npar
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
545 ioff2=ioff2+npar
546 END DO
547 END DO
548
549END SUBROUTINE qlsmq
550
551
563SUBROUTINE qlssq(aprod,A,s,roff,t)
564 USE mpqldec
565 USE mpdalc
566
567 ! cost[dot ops] ~= N*N*Nhr
568
569 IMPLICIT NONE
570 INTEGER(mpi) :: i
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
579 INTEGER(mpi) :: j
580 INTEGER(mpi) :: k
581 INTEGER(mpi) :: l
582 INTEGER(mpi) :: kn
583 INTEGER(mpl) :: length
584 INTEGER(mpi) :: nconb
585 INTEGER(mpi) :: nparb
586 REAL(mpd) :: vtAv
587 REAL(mpd), DIMENSION(:), ALLOCATABLE :: Av
588
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
593
594 INTERFACE
595 SUBROUTINE aprod(n,l,x,is,ie,y) ! y=A*x
596 USE mpdef
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)
603 END SUBROUTINE aprod
604 END INTERFACE
605 !$POMP INST BEGIN(qlssq)
606
607 length=npar
608 CALL mpalloc(av,length,'qlssq: A*v')
609
610 ioffp=0 ! parameter offset for block
611 DO ibpar=1,nblock ! parameter block
612 icoff=ioffblock(ibpar) ! constraint offset in parameter block
613 iclast=ioffblock(ibpar+1) ! last constraint in parameter block
614 nconb=iclast-icoff ! number of constraints in block
615 nparb=nparblock(ibpar) ! number of parameters in block
616 DO j=1,nconb
617 k=j
618 ! monitoring
619 IF(monpg>0) CALL monpgs(icoff+k)
620 IF (t) k=nconb+1-j
621 kn=nparb+k-nconb
622 ! expand row 'l' of matV into vecN
623 l=k+icoff
624 ! non-zero range (excluding 'kn')
625 ifirst=irangeparnz(1,l)
626 ilast=irangeparnz(2,l)
627 vecn(1:nparb)=0._mpd
628 vecn(ifirst:ilast)=matv(ioffrow(l)+1:ioffrow(l)+1+ilast-ifirst)
629 vecn(kn)=vecvk(k)
630 ! A*v
631 av(1:nparb)=0._mpd
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))
634 ! transformation
635 ! diagonal block
636 ! v^t*A*v
637 vtav=dot_product(vecn(ifirst:ilast),av(ifirst:ilast))+vecn(kn)*av(kn)
638 ! update
639 ! parallelize row loop
640 ! slot of 8 'I' for next idle thread
641 !$OMP PARALLEL DO &
642 !$OMP PRIVATE(IOFF2,ILASTI) &
643 !$OMP SCHEDULE(DYNAMIC,8)
644 DO i=1,kn
645 ioff2=roff(i+ioffp)+ioffp
646 ilasti=min(ilast,i)
647 ! correct with 2*(2v*vtAv*v^t - Av*v^t)
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))
650 END DO
651 !$OMP END PARALLEL DO
652
653 ! parallelize row loop
654 ! slot of 8 'I' for next idle thread
655 !$OMP PARALLEL DO &
656 !$OMP PRIVATE(IOFF2) &
657 !$OMP SCHEDULE(DYNAMIC,8)
658 DO i=ifirst,ilast
659 ioff2=roff(i+ioffp)+ioffp
660 ! correct with -2(Av*v^t)^t)
661 a(ioff2+1:ioff2+i)=a(ioff2+1:ioff2+i)-2.0_mpd*av(1:i)*vecn(i)
662 END DO
663 !$OMP END PARALLEL DO
664 ! i=kn, add secondary diagonal element
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))
667 ! off diagonal block
668 DO i=kn+1,nparb
669 ioff2=roff(i+ioffp)+ioffp
670 ! correct with -2Av*v^t
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)
673 END DO
674 END DO
675 ! update parameter offset
676 ioffp=ioffp+nparb
677 END DO
678
679 CALL mpdealloc(av)
680 !$POMP INST END(qlssq)
681
682END SUBROUTINE qlssq
683
695SUBROUTINE qlpssq(aprod,B,m,t)
696 USE mpqldec
697 USE mpdalc
698
699 ! cost[dot ops] ~= N*N*Nhr
700
701 IMPLICIT NONE
702 INTEGER(mpi) :: 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)
710 INTEGER(mpi) :: j
711 INTEGER(mpi) :: j2
712 INTEGER(mpi) :: k
713 INTEGER(mpi) :: k2
714 INTEGER(mpi) :: kn
715 INTEGER(mpi) :: kn2
716 INTEGER(mpi) :: l
717 INTEGER(mpi) :: l1
718 INTEGER(mpi) :: l2
719 INTEGER(mpl) :: length
720 INTEGER(mpi) :: mbnd
721 REAL(mpd) :: v2kn
722 REAL(mpd) :: vtAv
723 REAL(mpd) :: vtAvp
724 REAL(mpd) :: vtvp
725 REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecAv ! A*v
726 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matvtvp ! v^t*v'
727 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matvtAvp ! v^t*(A*v')
728 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matCoeff ! coefficients (d(A*v)=sum(c_i*v_i))
729 INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: irangeCoeff
730
731 INTEGER(mpi), INTENT(IN) :: m
732 REAL(mpd), INTENT(IN OUT) :: B(npar*m-(m*m-m)/2)
733 LOGICAL, INTENT(IN) :: t
734
735 INTERFACE
736 SUBROUTINE aprod(n,l,x,is,ie,y) ! y=A*x
737 USE mpdef
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)
744 END SUBROUTINE aprod
745 END INTERFACE
746 !$POMP INST BEGIN(qlpssq)
747
748 length=npar
749 CALL mpalloc(vecav,length,'qlpssq: A*v')
750 length=int(ncon,mpl)*int(ncon,mpl)
751 CALL mpalloc(matvtvp,length,"qlpssq: v^t*v'")
752 matvtvp=0._mpd
753 CALL mpalloc(matvtavp,length,"qlpssq: v^t*(A*v')")
754 matvtavp=0._mpd
755 CALL mpalloc(matcoeff,length,'qlpssq: coefficients')
756 matcoeff=0._mpd
757 length=ncon
758 CALL mpalloc(irangecoeff,2_mpl,length,'qlpssq: non zero coefficient range')
759
760 mbnd=max(0,m-1) ! band width without diagonal
761
762 DO j=1,ncon
763 k=j
764 ! monitoring
765 IF(monpg>0) CALL monpgs(k)
766 IF (t) k=ncon+1-j
767 kn=npar+k-ncon
768 ioff1=int(k-1,mpl)*int(ncon,mpl)
769 irangecoeff(1,k)=ncon
770 irangecoeff(2,k)=1
771 ! expand row 'k' of matV into vecN
772 ! non-zero range (excluding 'kn')
773 ifirst=irangeparnz(1,k)
774 ilast=irangeparnz(2,k)
775 vecn=0._mpd
776 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
777 vecn(kn)=vecvk(k)
778 ! transformation A*v
779 vecav(1:npar)=0._mpd
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))
782 ! products v^t*v'
783 DO j2=j+1,ncon
784 k2=j2
785 IF (t) k2=ncon+1-j2
786 kn2=npar+k2-ncon
787 ioff2=int(k2-1,mpl)*int(ncon,mpl)
788 ! non-zero range (excluding 'kn')
789 ifirst2=irangeparnz(1,k2)
790 ilast2=irangeparnz(2,k2)
791 v2kn=0._mpd
792 IF (kn >= ifirst2.AND.kn <= ilast2) v2kn=matv(ioffrow(k2)+1+kn-ifirst2)
793 ! overlap regions
794 l1=max(ifirst,ifirst2)
795 l2=min(ilast,ilast2)
796 vtvp=vecn(kn2)*vecvk(k2)+vecn(kn)*v2kn ! v^t*v'
797 IF (l1 <= l2) vtvp=vtvp+dot_product(vecn(l1:l2), &
798 matv(ioffrow(k2)+1+l1-ifirst2:ioffrow(k2)+1+l2-ifirst2))
799 ! significant term?
800 IF (abs(vtvp) > 16.0_mpd*epsilon(vtvp)) THEN
801 matvtvp(ioff1+k2)=vtvp
802 matvtvp(ioff2+k)=vtvp
803 END IF
804 END DO
805 matvtvp(ioff1+k)=1.0_mpd
806 ! products v^t*(A*v')
807 DO j2=1,j
808 k2=j2
809 IF (t) k2=ncon+1-j2
810 kn2=npar+k2-ncon
811 ! non-zero range (excluding 'kn')
812 ifirst2=irangeparnz(1,k2)
813 ilast2=irangeparnz(2,k2)
814 ! non-zero regions
815 matvtavp(ioff1+k2)=vecvk(k2)*vecav(kn2)+dot_product(vecav(ifirst2:ilast2), &
816 matv(ioffrow(k2)+1:ioffrow(k2)+1+ilast2-ifirst2)) ! v'^t*(A*v)
817 END DO
818 ! update with (initial) A*v
819 ioff2=0
820 vtav=matvtavp(ioff1+k)
821 DO i=1,kn
822 ! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
823 DO l=max(1,i-mbnd),i
824 ioff2=ioff2+1
825 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*vecn(i)*vtav-vecav(i))*vecn(l)-vecav(l)*vecn(i))
826 END DO
827 END DO
828 ! off diagonal block
829 DO i=kn+1,npar
830 ! correct with -2Av*v^t
831 DO l=max(1,i-mbnd),i
832 ioff2=ioff2+1
833 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*vecn(l)
834 END DO
835 END DO
836 END DO
837
838 ! corrections for A*v (as linear combination of v's)
839 DO j=1,ncon
840 k=j
841 IF (t) k=ncon+1-j
842 kn=npar+k-ncon
843 ioff1=int(k-1,mpl)*int(ncon,mpl)
844 ! expand row 'k' of matV into vecN
845 ! non-zero range (excluding 'kn')
846 ifirst=irangeparnz(1,k)
847 ilast=irangeparnz(2,k)
848 vecn=0._mpd
849 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
850 vecn(kn)=vecvk(k)
851 ! transformation (diagonal block)
852 l1=irangecoeff(1,k)
853 l2=irangecoeff(2,k)
854 ! diagonal block
855 ! v^t*A*v
856 vtav=matvtavp(ioff1+k)+dot_product(matcoeff(ioff1+l1:ioff1+l2),matvtvp(ioff1+l1:ioff1+l2))
857 ! expand correction to initial A*v
858 vecav(1:npar)=0._mpd
859 istat=0
860 DO k2=l1,l2
861 IF (matcoeff(ioff1+k2) == 0._mpd) cycle
862 if (istat(1)==0) istat(1)=k2
863 istat(2)=k2
864 istat(3)=istat(3)+1
865 kn2=npar+k2-ncon
866 ! expand row 'k2' of matV directly into vecAv
867 ! non-zero range (excluding 'kn')
868 ifirst2=irangeparnz(1,k2)
869 ilast2=irangeparnz(2,k2)
870 vecav(ifirst2:ilast2)=vecav(ifirst2:ilast2)+matcoeff(ioff1+k2)* &
871 matv(ioffrow(k2)+1:ioffrow(k2)+1+ilast2-ifirst2)
872 vecav(kn2)=vecav(kn2)+matcoeff(ioff1+k2)*vecvk(k2)
873 END DO
874 ! update
875 ioff2=0
876 DO i=1,kn
877 ! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
878 DO l=max(1,i-mbnd),i
879 ioff2=ioff2+1
880 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*vecn(i)*vtav-vecav(i))*vecn(l)-vecav(l)*vecn(i))
881 END DO
882 END DO
883 ! off diagonal block
884 DO i=kn+1,npar
885 ! correct with -2Av*v^t
886 DO l=max(1,i-mbnd),i
887 ioff2=ioff2+1
888 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*vecn(l)
889 END DO
890 END DO
891 ! correct A*v for the remainung v
892 DO j2=j+1,ncon
893 k2=j2
894 IF (t) k2=ncon+1-j2
895 kn2=npar+k2-ncon
896 ioff2=int(k2-1,mpl)*int(ncon,mpl)
897 vtvp=matvtvp(ioff1+k2) ! v^t*v'
898 ! non-zero regions
899 l1=irangecoeff(1,k2)
900 l2=irangecoeff(2,k2)
901 vtavp=matvtavp(ioff2+k)
902 IF (l1 <= l2) vtavp=vtavp+dot_product(matcoeff(ioff2+l1:ioff2+l2),matvtvp(ioff1+l1:ioff1+l2)) ! v^t*(A*v')
903 l1=min(l1,k)
904 l2=max(l2,k)
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
910 END IF
911 irangecoeff(1,k2)=l1
912 irangecoeff(2,k2)=l2
913 END DO
914 END DO
915
916 CALL mpdealloc(irangecoeff)
917 CALL mpdealloc(matcoeff)
918 CALL mpdealloc(matvtavp)
919 CALL mpdealloc(matvtvp)
920 CALL mpdealloc(vecav)
921 !$POMP INST END(qlpssq)
922
923END SUBROUTINE qlpssq
924
925
933SUBROUTINE qlgete(emin,emax)
934 USE mpqldec
935
936 IMPLICIT NONE
937 INTEGER(mpi) :: i
938 INTEGER(mpi) :: ibpar
939 INTEGER(mpi) :: icoff
940 INTEGER(mpi) :: iclast
941 INTEGER(mpl) :: idiag
942
943 REAL(mpd), INTENT(OUT) :: emin
944 REAL(mpd), INTENT(OUT) :: emax
945
946 emax=matl(1)
947 emin=emax
948 DO ibpar=1,nblock ! parameter block
949 icoff=ioffblock(ibpar) ! constraint offset in parameter block
950 iclast=ioffblock(ibpar+1) ! last constraint in parameter block
951 idiag=int(ncon,mpl)*int(icoff,mpl)+1
952 DO i=icoff+1,iclast
953 IF (abs(emax) < abs(matl(idiag))) emax=matl(idiag)
954 IF (abs(emin) > abs(matl(idiag))) emin=matl(idiag)
955 idiag=idiag+ncon+1
956 END DO
957 END DO
958
959END SUBROUTINE qlgete
960
961
969SUBROUTINE qlbsub(d,y)
970 USE mpqldec
971
972 IMPLICIT NONE
973 INTEGER(mpi) :: icoff
974 INTEGER(mpi) :: iclast
975 INTEGER(mpl) :: idiag
976 INTEGER(mpi) :: k
977 INTEGER(mpi) :: nconb
978
979 REAL(mpd), INTENT(IN) :: d(ncon)
980 REAL(mpd), INTENT(OUT) :: y(ncon)
981
982 ! solve L*y=d by forward substitution
983 icoff=ioffblock(iblock) ! constraint offset in parameter block
984 iclast=ioffblock(iblock+1) ! last constraint in parameter block
985 nconb=iclast-icoff ! number of constraints in block
986 idiag=int(ncon,mpl)*int(iclast-1,mpl)+nconb
987 DO k=nconb,1,-1
988 y(k)=(d(k)-dot_product(matl(idiag+1:idiag+nconb-k),y(k+1:nconb)))/matl(idiag)
989 idiag=idiag-ncon-1
990 END DO
991
992END SUBROUTINE qlbsub
993
996SUBROUTINE qlsetb(ib)
997 USE mpqldec
998
999 IMPLICIT NONE
1000 INTEGER(mpi), INTENT(IN) :: ib
1001
1002 iblock=ib
1003
1004END SUBROUTINE qlsetb
1005
1008SUBROUTINE qldump()
1009 USE mpqldec
1010
1011 IMPLICIT NONE
1012 INTEGER(mpi) :: i
1013 INTEGER(mpi) :: ifirst
1014 INTEGER(mpi) :: ilast
1015 INTEGER(mpl) :: ioff1
1016 INTEGER(mpl) :: ioff2
1017 INTEGER(mpi) :: istat(6)
1018 INTEGER(mpi) :: j
1019 INTEGER(mpi) :: kn
1020 REAL(mpd) :: v1
1021 REAL(mpd) :: v2
1022 REAL(mpd) :: v3
1023 REAL(mpd) :: v4
1024
1025 print *
1026 ioff1=0
1027 ioff2=0
1028
1029 DO i=1, ncon
1030 kn=npar-ncon+i
1031 istat=0
1032 v1=0.;v2=0.;v3=0.;v4=0.
1033 ! expand row 'i' of matV into vecN
1034 ! non-zero range (excluding 'kn')
1035 ifirst=irangeparnz(1,i)
1036 ilast=irangeparnz(2,i)
1037 vecn=0._mpd
1038 vecn(ifirst:ilast)=matv(ioffrow(i)+1:ioffrow(i)+1+ilast-ifirst)
1039 DO j=1,npar+i-ncon
1040 IF (vecn(j) /= 0.0_mpd) THEN
1041 v2=vecn(j)
1042 IF (istat(3) == 0) THEN
1043 istat(1)=j
1044 v1=v2
1045 END IF
1046 istat(2)=j
1047 istat(3)=istat(3)+1
1048 END IF
1049 END DO
1050 ioff1=ioff1+npar
1051 DO j=1,ncon
1052 IF (matl(ioff2+j) /= 0.0_mpd) THEN
1053 v4=matl(ioff2+j)
1054 IF (istat(6) == 0) THEN
1055 istat(4)=j
1056 v3=v4
1057 END IF
1058 istat(5)=j
1059 istat(6)=istat(6)+1
1060 END IF
1061 END DO
1062 ioff2=ioff2+ncon
1063 print 100, i, istat, v1, v2, v3, v4, vecvk(i), irangeparnz(:,i)
1064 END DO
1065 print *
1066100 FORMAT(" qldump",7i8,5g13.5,2i8)
1067
1068END SUBROUTINE qldump
allocate array
Definition: mpdalc.f90:36
deallocate array
Definition: mpdalc.f90:42
subroutine monpgs(i)
Progress monitoring.
Definition: mpmon.f90:68
subroutine qlmrq(x, m, t)
Multiply right by Q(t).
Definition: mpqldec.f90:456
subroutine qldump()
Print statistics.
Definition: mpqldec.f90:1009
subroutine qlsmq(x, t)
Similarity transformation by Q(t).
Definition: mpqldec.f90:504
subroutine qlpssq(aprod, B, m, t)
Partial similarity transformation by Q(t).
Definition: mpqldec.f90:696
subroutine qldecb(a, bpar, bcon, rcon)
QL decomposition (for disjoint block matrix).
Definition: mpqldec.f90:216
subroutine qldec(a)
QL decomposition (as single block).
Definition: mpqldec.f90:122
subroutine qlmlq(x, m, t)
Multiply left by Q(t) (per block).
Definition: mpqldec.f90:395
subroutine qlsetb(ib)
Set block.
Definition: mpqldec.f90:997
subroutine qlbsub(d, y)
Backward substitution (per block).
Definition: mpqldec.f90:970
subroutine qlini(n, m, l, s, k)
Initialize QL decomposition.
Definition: mpqldec.f90:58
subroutine qlgete(emin, emax)
Get eigenvalues.
Definition: mpqldec.f90:934
subroutine qlssq(aprod, A, s, roff, t)
Similarity transformation by Q(t).
Definition: mpqldec.f90:564
(De)Allocate vectors and arrays.
Definition: mpdalc.f90:24
Definition of constants.
Definition: mpdef.f90:24
integer, parameter mpd
double precision
Definition: mpdef.f90:38
QL data.
Definition: mpqldec.f90:27
integer(mpi) ncon
number of constraints
Definition: mpqldec.f90:32
integer(mpi) iblock
active block
Definition: mpqldec.f90:35
real(mpd), dimension(:), allocatable vecvk
secondary diagonal of matV (including last element)
Definition: mpqldec.f90:38
integer(mpi), dimension(:,:), allocatable irangeparnz
range for non zero part (except vecVk)
Definition: mpqldec.f90:45
integer(mpi) monpg
flag for progress monitoring
Definition: mpqldec.f90:36
integer(mpl), dimension(:), allocatable ioffrow
row offsets in matV (for constrint block)
Definition: mpqldec.f90:43
integer(mpl) matsize
size of contraints matrix
Definition: mpqldec.f90:34
real(mpd), dimension(:), allocatable matv
unit normals (v_i) of Householder reflectors
Definition: mpqldec.f90:37
real(mpd), dimension(:), allocatable matl
lower diagonal matrix L
Definition: mpqldec.f90:39
integer(mpi), dimension(:), allocatable ioffpar
parameter number offsets for matV ( " )
Definition: mpqldec.f90:44
integer(mpi) nblock
number of blocks
Definition: mpqldec.f90:33
real(mpd), dimension(:), allocatable vecn
normal vector
Definition: mpqldec.f90:40
integer(mpi) npar
number of parameters
Definition: mpqldec.f90:31
integer(mpi), dimension(:), allocatable ioffblock
block offset (1.
Definition: mpqldec.f90:42
integer(mpi), dimension(:), allocatable nparblock
number of parameters in block
Definition: mpqldec.f90:41