Millepede-II V04-17-04
mpqldec.f90
Go to the documentation of this file.
1
25
26#ifdef SCOREP_USER_ENABLE
27#include "scorep/SCOREP_User.inc"
28#endif
29
31MODULE mpqldec
32 USE mpdef
33 IMPLICIT NONE
34
35 INTEGER(mpi) :: npar
36 INTEGER(mpi) :: ncon
37 INTEGER(mpi) :: nblock
38 INTEGER(mpl) :: matsize
39 INTEGER(mpi) :: iblock
40 INTEGER(mpi) :: monpg
41 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matv
42 REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecvk
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
49 INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: irangeparnz
50
51END MODULE mpqldec
52
61SUBROUTINE qlini(n,m,l,s,k)
62 USE mpqldec
63 USE mpdalc
64
65 IMPLICIT NONE
66 INTEGER(mpl) :: length
67
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
73
74 npar=n
75 ncon=m
76 nblock=l
77 matsize=s
78 iblock=1
79 monpg=k
80 ! allocate
81 length=matsize
82 !print *, ' full length (V)', length
83 CALL mpalloc(matv,length,'QLDEC: V')
84 matv=0.
85 length=int(ncon,mpl)*int(ncon,mpl)
86 CALL mpalloc(matl,length,'QLDEC: L')
87 matl=0.
88 length=npar
89 CALL mpalloc(vecn,length,'QLDEC: v')
90 length=ncon
91 CALL mpalloc(vecvk,length,'QLDEC: sec. diag(V)')
92 vecvk=0.
93 CALL mpalloc(ioffpar,length,'QLDEC: parameter offsets (V)')
94 ioffpar=0
95 CALL mpalloc(irangeparnz,2_mpl,length,'QLDEC: parameter non zero range (V)')
96 length=ncon+1
97 CALL mpalloc(ioffrow,length,'QLDEC: row offsets (V)')
98 ioffrow=0
99 length=nblock
100 CALL mpalloc(nparblock,length,'QLDEC: npar in block')
101 nparblock=0
102 length=nblock+1
103 CALL mpalloc(ioffblock,length,'QLDEC: ioff for block')
104 ioffblock=0
105END SUBROUTINE qlini
106
107! 141217 C. Kleinwort, DESY-FH1
125SUBROUTINE qldec(a)
126 USE mpqldec
127 USE mpdalc
128
129 ! cost[dot ops] ~= Npar*Ncon*Ncon
130
131 IMPLICIT NONE
132 INTEGER(mpi) :: i
133 INTEGER(mpl) :: ioff1
134 INTEGER(mpl) :: ioff2
135 INTEGER(mpl) :: ioff3
136 INTEGER(mpi) :: k
137 INTEGER(mpi) :: kn
138 INTEGER(mpl) :: length
139 REAL(mpd) :: nrm
140 REAL(mpd) :: sp
141
142 REAL(mpd), INTENT(IN) :: a(matsize)
143
144 ! prepare
145 vecvk=0._mpd
146 length=int(npar,mpl)*int(ncon,mpl)
147 matv=a(1:length)
148 matl=0.0_mpd
149 ! implemented as single block
150 nblock=1
151 nparblock(1)=npar
152 ioffblock(2)=ncon
153 DO k=1,ncon
154 ioffrow(k+1)=ioffrow(k)+npar
155 END DO
156
157 ! Householder procedure
158 DO k=ncon,1,-1
159 ! monitoring
160 IF(monpg>0) CALL monpgs(ncon+1-k)
161 kn=npar+k-ncon
162 ! column offset
163 ioff1=int(k-1,mpl)*int(npar,mpl)
164 ! get column
165 vecn(1:kn)=matv(ioff1+1:ioff1+kn)
166 nrm = sqrt(dot_product(vecn(1:kn),vecn(1:kn)))
167 IF (nrm == 0.0_mpd) cycle
168 !
169 IF (vecn(kn) >= 0.0_mpd) THEN
170 vecn(kn)=vecn(kn)+nrm
171 ELSE
172 vecn(kn)=vecn(kn)-nrm
173 END IF
174 ! create normal vector
175 nrm = sqrt(dot_product(vecn(1:kn),vecn(1:kn)))
176 vecn(1:kn)=vecn(1:kn)/nrm
177 ! transformation
178 ioff2=0
179 DO i=1,k
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
182 ioff2=ioff2+npar
183 END DO
184 ! store column of L
185 ioff3=int(k-1,mpl)*int(ncon,mpl)
186 matl(ioff3+k:ioff3+ncon)=matv(ioff1+kn:ioff1+npar)
187 ! store normal vector
188 matv(ioff1+1:ioff1+kn-1)=vecn(1:kn-1)
189 matv(ioff1+kn:ioff1+npar)=0.0_mpd
190 irangeparnz(1,k)=1
191 irangeparnz(2,k)=kn-1
192 ! store secondary diagonal
193 vecvk(k)=vecn(kn)
194 END DO
195
196END SUBROUTINE qldec
197
198! 190312 C. Kleinwort, DESY-BELLE
219SUBROUTINE qldecb(a,bpar,bcon,rcon)
220 USE mpqldec
221 USE mpdalc
222
223 ! cost[dot ops] ~= Npar*Ncon*Ncon
224
225 IMPLICIT NONE
226 INTEGER(mpi) :: i
227 INTEGER(mpi) :: ibcon
228 INTEGER(mpi) :: iblast
229 INTEGER(mpi) :: iboff
230 INTEGER(mpi) :: ibpar
231 INTEGER(mpi) :: ifirst
232 INTEGER(mpi) :: ilast
233 INTEGER(mpi) :: in
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
241 INTEGER(mpi) :: k
242 INTEGER(mpi) :: k1
243 INTEGER(mpi) :: kn
244 INTEGER(mpi) :: ncb
245 INTEGER(mpi) :: ncol
246 INTEGER(mpi) :: npb
247 REAL(mpd) :: nrm
248 REAL(mpd) :: sp
249
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)
254
255 !$POMP INST BEGIN(qldecb)
256#ifdef SCOREP_USER_ENABLE
257 scorep_user_region_by_name_begin("UR_qldecb", scorep_user_region_type_common)
258#endif
259 ! prepare
260 vecvk=0.0_mpd
261 matv=a(1:matsize)
262 matl=0.0_mpd
263
264 ! prepare offsets
265 icoff=0
266 DO ibpar=1,nblock ! parameter block
267 iclast=icoff
268 DO ibcon=bpar(2,ibpar)+1, bpar(2,ibpar+1)! constraint block
269 ncb=bcon(1,ibcon+1)-bcon(1,ibcon) ! number of constraints in constraint block
270 npb=bcon(3,ibcon)+1-bcon(2,ibcon) ! number of parameters in constraint block
271 ifirst=bcon(2,ibcon)
272 ilast=bcon(3,ibcon)
273 DO i=bcon(1,ibcon),bcon(1,ibcon+1)-1
274 ! non-zero range: first, last parameter
275 irangeparnz(1,i)=rcon(1,i)
276 irangeparnz(2,i)=rcon(2,i)
277 ! storage: parameter, row offset
278 ioffpar(i)=rcon(3,i)-1
279 ioffrow(i+1)=ioffrow(i)+rcon(4,i)-ioffpar(i)
280 END DO
281 iclast=iclast+ncb
282 END DO
283 ! set up matL
284 iblast=bpar(1,ibpar+1) ! last parameter in parameter block
285 DO k=icoff+1,iclast
286 kn=iblast+k-iclast
287 ioff1=ioffrow(k)
288 npb=int(ioffrow(k+1)-ioff1,mpi)
289 ! size of overlap
290 ncol=ioffpar(k)+npb-kn
291 IF (ncol >= 0) THEN
292 ioff3=int(k-1,mpl)*int(ncon,mpl)
293 matl(ioff3+iclast-ncol-icoff:ioff3+iclast-icoff)=matv(ioff1+npb-ncol:ioff1+npb)
294 END IF
295 END DO
296 icoff=iclast
297 nparblock(ibpar)=bpar(1,ibpar+1)-bpar(1,ibpar)
298 ioffblock(ibpar+1)=icoff
299 END DO
300
301 DO ibpar=1,nblock ! parameter block
302 iboff=bpar(1,ibpar) ! last offset in parameter block
303 iblast=bpar(1,ibpar+1) ! last parameter in parameter block
304 icoff=ioffblock(ibpar) ! constraint offset in parameter block
305 iclast=ioffblock(ibpar+1) ! last constraint in parameter block
306 IF(iclast <= icoff) cycle ! no constraints
307 ibcon=bpar(2,ibpar+1) ! start with last constraint block
308 k1=bcon(1,ibcon) ! first constraint in block
309 ! Householder procedure
310 DO k=iclast,icoff+1,-1
311 ! monitoring
312 IF(monpg>0) CALL monpgs(ncon+1-k)
313 kn=iblast+k-iclast
314 ! different constraint block?
315 IF (k < k1) THEN
316 ibcon=ibcon-1
317 k1=bcon(1,ibcon)
318 END IF
319 ! parameter offset
320 ipoff=ioffpar(k)
321 ! index if first non-zero parameter
322 ifirst=ipoff+1
323 IF (ifirst > kn) cycle
324 ! column offset
325 ioff1=ioffrow(k)
326 ! number of parameter
327 npb=int(ioffrow(k+1)-ioff1,mpi)
328 ! index of last parameter
329 iplast=ioffpar(k)+npb
330 ! index of last used parameter
331 ilast=min(iplast,kn)
332 ! number of used columns
333 ncol=ilast-ipoff
334 ! get column
335 in=iblast+k1-iclast
336 vecn(in:kn)=0.0_mpd
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
340 !
341 IF (vecn(kn) >= 0.0_mpd) THEN
342 vecn(kn)=vecn(kn)+nrm
343 ELSE
344 vecn(kn)=vecn(kn)-nrm
345 END IF
346 ! create normal vector
347 IF (ilast < kn) THEN
348 nrm = sqrt(dot_product(vecn(ifirst:ilast),vecn(ifirst:ilast))+vecn(kn)*vecn(kn))
349 vecn(ifirst:ilast)=vecn(ifirst:ilast)/nrm
350 vecn(kn)=vecn(kn)/nrm
351 ELSE
352 nrm = sqrt(dot_product(vecn(ifirst:ilast),vecn(ifirst:ilast)))
353 vecn(ifirst:ilast)=vecn(ifirst:ilast)/nrm
354 END IF
355 ! update L too
356 ioff3=int(k1-2,mpl)*int(ncon,mpl)
357 ! transformation
358 DO i=k1,k
359 ioff3=ioff3+ncon
360 IF (irangeparnz(1,k) > irangeparnz(2,i)) cycle ! no overlap
361 ioff2=ioffrow(i)+ioffpar(k)-ioffpar(i)
362 sp=dot_product(vecn(ifirst:ilast),matv(ioff2+1:ioff2+ncol))
363 IF (sp /= 0.0_mpd) THEN
364 ! update matV
365 matv(ioff2+1:ioff2+ncol)=matv(ioff2+1:ioff2+ncol)-2.0_mpd*vecn(ifirst:ilast)*sp
366 ! update matL
367 in=iblast+i-iclast
368 matl(ioff3+i-icoff:ioff3+k-icoff)=matl(ioff3+i-icoff:ioff3+k-icoff)-2.0_mpd*vecn(in:kn)*sp
369 ! update non zero range
370 irangeparnz(1,i)=min(irangeparnz(1,i),irangeparnz(1,k))
371 irangeparnz(2,i)=max(irangeparnz(2,i),irangeparnz(2,k))
372 END IF
373 END DO
374 ! store secondary diagonal
375 vecvk(icoff+k)=vecn(kn)
376 ! store normal vector (non zero part)
377 ifirst=irangeparnz(1,k)
378 ilast=min(irangeparnz(2,k),kn-1)
379 ncol=ilast-ifirst+1
380 matv(ioff1+1:ioff1+ncol)=vecn(ifirst:ilast)
381 matv(ioff1+ncol+1:ioff1+npb)=0.0_mpd
382 ! local to parameter block
383 irangeparnz(1,k)=ifirst-iboff
384 irangeparnz(2,k)=ilast-iboff
385 END DO
386 END DO
387#ifdef SCOREP_USER_ENABLE
388 scorep_user_region_by_name_end("UR_qldecb")
389#endif
390 !$POMP INST END(qldecb)
391
392
393END SUBROUTINE qldecb
394
395
404SUBROUTINE qlmlq(x,m,t)
405 USE mpqldec
406
407 ! cost[dot ops] ~= N*M*Nhr
408
409 IMPLICIT NONE
410 INTEGER(mpi) :: i
411 INTEGER(mpi) :: icoff
412 INTEGER(mpi) :: iclast
413 INTEGER(mpi) :: ifirst
414 INTEGER(mpi) :: ilast
415 INTEGER(mpl) :: ioff2
416 INTEGER(mpi) :: j
417 INTEGER(mpi) :: k
418 INTEGER(mpi) :: l
419 INTEGER(mpi) :: kn
420 INTEGER(mpi) :: nconb
421 INTEGER(mpi) :: nparb
422 REAL(mpd) :: sp
423
424 INTEGER(mpi), INTENT(IN) :: m
425 REAL(mpd), INTENT(IN OUT) :: x(INT(npar,mpl)*INT(m,mpl))
426 LOGICAL, INTENT(IN) :: t
427
428 icoff=ioffblock(iblock) ! constraint offset in parameter block
429 iclast=ioffblock(iblock+1) ! last constraint in parameter block
430 nconb=iclast-icoff ! number of constraints in block
431 nparb=nparblock(iblock) ! number of parameters in block
432 DO j=1,nconb
433 k=j
434 IF (t) k=nconb+1-j
435 kn=nparb+k-nconb
436 ! expand row 'l' of matV into vecN
437 l=k+icoff
438 ! non-zero range (excluding 'kn')
439 ifirst=irangeparnz(1,l)
440 ilast=irangeparnz(2,l)
441 vecn(1:nparb)=0._mpd
442 vecn(ifirst:ilast)=matv(ioffrow(l)+1:ioffrow(l)+1+ilast-ifirst)
443 vecn(kn)=vecvk(l)
444 ! transformation
445 ioff2=0
446 DO i=1,m
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
450 ioff2=ioff2+nparb
451 END DO
452 END DO
453
454END SUBROUTINE qlmlq
455
456
465SUBROUTINE qlmrq(x,m,t)
466 USE mpqldec
467
468 ! cost[dot ops] ~= N*M*Nhr
469
470 IMPLICIT NONE
471 INTEGER(mpi) :: i
472 INTEGER(mpl) :: iend
473 INTEGER(mpi) :: ifirst
474 INTEGER(mpi) :: ilast
475 INTEGER(mpi) :: j
476 INTEGER(mpi) :: k
477 INTEGER(mpi) :: kn
478 REAL(mpd) :: sp
479
480 INTEGER(mpi), INTENT(IN) :: m
481 REAL(mpd), INTENT(IN OUT) :: x(INT(m,mpl)*INT(npar,mpl))
482 LOGICAL, INTENT(IN) :: t
483
484 DO j=1,ncon
485 k=j
486 IF (.not.t) k=ncon+1-j
487 kn=npar+k-ncon
488 ! expand row 'k' of matV into vecN
489 ! non-zero range (excluding 'kn')
490 ifirst=irangeparnz(1,k)
491 ilast=irangeparnz(2,k)
492 vecn=0._mpd
493 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
494 vecn(kn)=vecvk(k)
495 ! transformation
496 iend=m*kn
497 DO i=1,npar
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
500 END DO
501 END DO
502
503END SUBROUTINE qlmrq
504
505
513SUBROUTINE qlsmq(x,t)
514 USE mpqldec
515
516 ! cost[dot ops] ~= N*N*Nhr
517
518 IMPLICIT NONE
519 INTEGER(mpi) :: i
520 INTEGER(mpl) :: ioff2
521 INTEGER(mpl) :: iend
522 INTEGER(mpi) :: ifirst
523 INTEGER(mpi) :: ilast
524 INTEGER(mpi) :: j
525 INTEGER(mpi) :: k
526 INTEGER(mpi) :: kn
527 REAL(mpd) :: sp
528
529 REAL(mpd), INTENT(IN OUT) :: x(INT(npar,mpl)*INT(npar,mpl))
530 LOGICAL, INTENT(IN) :: t
531
532 DO j=1,ncon
533 ! monitoring
534 IF(monpg>0) CALL monpgs(j)
535 k=j
536 IF (t) k=ncon+1-j
537 kn=npar+k-ncon
538 ! expand row 'k' of matV into vecN
539 ! non-zero range (excluding 'kn')
540 ifirst=irangeparnz(1,k)
541 ilast=irangeparnz(2,k)
542 vecn=0._mpd
543 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
544 vecn(kn)=vecvk(k)
545 ! transformation
546 iend=int(npar,mpl)*int(kn,mpl)
547 DO i=1,npar
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
550 END DO
551 ioff2=0
552 DO i=1,npar
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
555 ioff2=ioff2+npar
556 END DO
557 END DO
558
559END SUBROUTINE qlsmq
560
561
573SUBROUTINE qlssq(aprod,A,s,roff,t)
574 USE mpqldec
575 USE mpdalc
576
577 ! cost[dot ops] ~= N*N*Nhr
578
579 IMPLICIT NONE
580 INTEGER(mpi) :: i
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
589 INTEGER(mpi) :: j
590 INTEGER(mpi) :: k
591 INTEGER(mpi) :: l
592 INTEGER(mpi) :: kn
593 INTEGER(mpl) :: length
594 INTEGER(mpi) :: nconb
595 INTEGER(mpi) :: nparb
596 REAL(mpd) :: vtAv
597 REAL(mpd), DIMENSION(:), ALLOCATABLE :: Av
598
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
603
604 INTERFACE
605 SUBROUTINE aprod(n,l,x,is,ie,y) ! y=A*x
606 USE mpdef
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)
613 END SUBROUTINE aprod
614 END INTERFACE
615 !$POMP INST BEGIN(qlssq)
616#ifdef SCOREP_USER_ENABLE
617 scorep_user_region_by_name_begin("UR_qlssq", scorep_user_region_type_common)
618#endif
619
620 length=npar
621 CALL mpalloc(av,length,'qlssq: A*v')
622
623 ioffp=0 ! parameter offset for block
624 DO ibpar=1,nblock ! parameter block
625 icoff=ioffblock(ibpar) ! constraint offset in parameter block
626 iclast=ioffblock(ibpar+1) ! last constraint in parameter block
627 nconb=iclast-icoff ! number of constraints in block
628 nparb=nparblock(ibpar) ! number of parameters in block
629 DO j=1,nconb
630 k=j
631 ! monitoring
632 IF(monpg>0) CALL monpgs(icoff+k)
633 IF (t) k=nconb+1-j
634 kn=nparb+k-nconb
635 ! expand row 'l' of matV into vecN
636 l=k+icoff
637 ! non-zero range (excluding 'kn')
638 ifirst=irangeparnz(1,l)
639 ilast=irangeparnz(2,l)
640 vecn(1:nparb)=0._mpd
641 vecn(ifirst:ilast)=matv(ioffrow(l)+1:ioffrow(l)+1+ilast-ifirst)
642 vecn(kn)=vecvk(k)
643 ! A*v
644 av(1:nparb)=0._mpd
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))
647 ! transformation
648 ! diagonal block
649 ! v^t*A*v
650 vtav=dot_product(vecn(ifirst:ilast),av(ifirst:ilast))+vecn(kn)*av(kn)
651 ! update
652 ! parallelize row loop
653 ! slot of 8 'I' for next idle thread
654 !$OMP PARALLEL DO &
655 !$OMP PRIVATE(IOFF2,ILASTI) &
656 !$OMP SCHEDULE(DYNAMIC,8)
657 DO i=1,kn
658 ioff2=roff(i+ioffp)+ioffp
659 ilasti=min(ilast,i)
660 ! correct with 2*(2v*vtAv*v^t - Av*v^t)
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))
663 END DO
664 !$OMP END PARALLEL DO
665
666 ! parallelize row loop
667 ! slot of 8 'I' for next idle thread
668 !$OMP PARALLEL DO &
669 !$OMP PRIVATE(IOFF2) &
670 !$OMP SCHEDULE(DYNAMIC,8)
671 DO i=ifirst,ilast
672 ioff2=roff(i+ioffp)+ioffp
673 ! correct with -2(Av*v^t)^t)
674 a(ioff2+1:ioff2+i)=a(ioff2+1:ioff2+i)-2.0_mpd*av(1:i)*vecn(i)
675 END DO
676 !$OMP END PARALLEL DO
677 ! i=kn, add secondary diagonal element
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))
680 ! off diagonal block
681 DO i=kn+1,nparb
682 ioff2=roff(i+ioffp)+ioffp
683 ! correct with -2Av*v^t
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)
686 END DO
687 END DO
688 ! update parameter offset
689 ioffp=ioffp+nparb
690 END DO
691
692 CALL mpdealloc(av)
693#ifdef SCOREP_USER_ENABLE
694 scorep_user_region_by_name_end("UR_qlssq")
695#endif
696 !$POMP INST END(qlssq)
697
698END SUBROUTINE qlssq
699
711SUBROUTINE qlpssq(aprod,B,m,t)
712 USE mpqldec
713 USE mpdalc
714
715 ! cost[dot ops] ~= N*N*Nhr
716
717 IMPLICIT NONE
718 INTEGER(mpi) :: i
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)
726 INTEGER(mpi) :: j
727 INTEGER(mpi) :: j2
728 INTEGER(mpi) :: k
729 INTEGER(mpi) :: k2
730 INTEGER(mpi) :: kn
731 INTEGER(mpi) :: kn2
732 INTEGER(mpi) :: l
733 INTEGER(mpi) :: l1
734 INTEGER(mpi) :: l2
735 INTEGER(mpl) :: length
736 INTEGER(mpi) :: mbnd
737 REAL(mpd) :: v2kn
738 REAL(mpd) :: vtAv
739 REAL(mpd) :: vtAvp
740 REAL(mpd) :: vtvp
741 REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecAv ! A*v
742 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matvtvp ! v^t*v'
743 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matvtAvp ! v^t*(A*v')
744 REAL(mpd), DIMENSION(:), ALLOCATABLE :: matCoeff ! coefficients (d(A*v)=sum(c_i*v_i))
745 INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: irangeCoeff
746
747 INTEGER(mpi), INTENT(IN) :: m
748 REAL(mpd), INTENT(IN OUT) :: B(npar*m-(m*m-m)/2)
749 LOGICAL, INTENT(IN) :: t
750
751 INTERFACE
752 SUBROUTINE aprod(n,l,x,is,ie,y) ! y=A*x
753 USE mpdef
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)
760 END SUBROUTINE aprod
761 END INTERFACE
762 !$POMP INST BEGIN(qlpssq)
763#ifdef SCOREP_USER_ENABLE
764 scorep_user_region_by_name_begin("UR_qlpssq", scorep_user_region_type_common)
765#endif
766
767 length=npar
768 CALL mpalloc(vecav,length,'qlpssq: A*v')
769 length=int(ncon,mpl)*int(ncon,mpl)
770 CALL mpalloc(matvtvp,length,"qlpssq: v^t*v'")
771 matvtvp=0._mpd
772 CALL mpalloc(matvtavp,length,"qlpssq: v^t*(A*v')")
773 matvtavp=0._mpd
774 CALL mpalloc(matcoeff,length,'qlpssq: coefficients')
775 matcoeff=0._mpd
776 length=ncon
777 CALL mpalloc(irangecoeff,2_mpl,length,'qlpssq: non zero coefficient range')
778
779 mbnd=max(0,m-1) ! band width without diagonal
780
781 DO j=1,ncon
782 k=j
783 ! monitoring
784 IF(monpg>0) CALL monpgs(k)
785 IF (t) k=ncon+1-j
786 kn=npar+k-ncon
787 ioff1=int(k-1,mpl)*int(ncon,mpl)
788 irangecoeff(1,k)=ncon
789 irangecoeff(2,k)=1
790 ! expand row 'k' of matV into vecN
791 ! non-zero range (excluding 'kn')
792 ifirst=irangeparnz(1,k)
793 ilast=irangeparnz(2,k)
794 vecn=0._mpd
795 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
796 vecn(kn)=vecvk(k)
797 ! transformation A*v
798 vecav(1:npar)=0._mpd
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))
801 ! products v^t*v'
802 DO j2=j+1,ncon
803 k2=j2
804 IF (t) k2=ncon+1-j2
805 kn2=npar+k2-ncon
806 ioff2=int(k2-1,mpl)*int(ncon,mpl)
807 ! non-zero range (excluding 'kn')
808 ifirst2=irangeparnz(1,k2)
809 ilast2=irangeparnz(2,k2)
810 v2kn=0._mpd
811 IF (kn >= ifirst2.AND.kn <= ilast2) v2kn=matv(ioffrow(k2)+1+kn-ifirst2)
812 ! overlap regions
813 l1=max(ifirst,ifirst2)
814 l2=min(ilast,ilast2)
815 vtvp=vecn(kn2)*vecvk(k2)+vecn(kn)*v2kn ! v^t*v'
816 IF (l1 <= l2) vtvp=vtvp+dot_product(vecn(l1:l2), &
817 matv(ioffrow(k2)+1+l1-ifirst2:ioffrow(k2)+1+l2-ifirst2))
818 ! significant term?
819 IF (abs(vtvp) > 16.0_mpd*epsilon(vtvp)) THEN
820 matvtvp(ioff1+k2)=vtvp
821 matvtvp(ioff2+k)=vtvp
822 END IF
823 END DO
824 matvtvp(ioff1+k)=1.0_mpd
825 ! products v^t*(A*v')
826 DO j2=1,j
827 k2=j2
828 IF (t) k2=ncon+1-j2
829 kn2=npar+k2-ncon
830 ! non-zero range (excluding 'kn')
831 ifirst2=irangeparnz(1,k2)
832 ilast2=irangeparnz(2,k2)
833 ! non-zero regions
834 matvtavp(ioff1+k2)=vecvk(k2)*vecav(kn2)+dot_product(vecav(ifirst2:ilast2), &
835 matv(ioffrow(k2)+1:ioffrow(k2)+1+ilast2-ifirst2)) ! v'^t*(A*v)
836 END DO
837 ! update with (initial) A*v
838 ioff2=0
839 vtav=matvtavp(ioff1+k)
840 DO i=1,kn
841 ! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
842 DO l=max(1,i-mbnd),i
843 ioff2=ioff2+1
844 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*vecn(i)*vtav-vecav(i))*vecn(l)-vecav(l)*vecn(i))
845 END DO
846 END DO
847 ! off diagonal block
848 DO i=kn+1,npar
849 ! correct with -2Av*v^t
850 DO l=max(1,i-mbnd),i
851 ioff2=ioff2+1
852 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*vecn(l)
853 END DO
854 END DO
855 END DO
856
857 ! corrections for A*v (as linear combination of v's)
858 DO j=1,ncon
859 k=j
860 IF (t) k=ncon+1-j
861 kn=npar+k-ncon
862 ioff1=int(k-1,mpl)*int(ncon,mpl)
863 ! expand row 'k' of matV into vecN
864 ! non-zero range (excluding 'kn')
865 ifirst=irangeparnz(1,k)
866 ilast=irangeparnz(2,k)
867 vecn=0._mpd
868 vecn(ifirst:ilast)=matv(ioffrow(k)+1:ioffrow(k)+1+ilast-ifirst)
869 vecn(kn)=vecvk(k)
870 ! transformation (diagonal block)
871 l1=irangecoeff(1,k)
872 l2=irangecoeff(2,k)
873 ! diagonal block
874 ! v^t*A*v
875 vtav=matvtavp(ioff1+k)+dot_product(matcoeff(ioff1+l1:ioff1+l2),matvtvp(ioff1+l1:ioff1+l2))
876 ! expand correction to initial A*v
877 vecav(1:npar)=0._mpd
878 istat=0
879 DO k2=l1,l2
880 IF (matcoeff(ioff1+k2) == 0._mpd) cycle
881 if (istat(1)==0) istat(1)=k2
882 istat(2)=k2
883 istat(3)=istat(3)+1
884 kn2=npar+k2-ncon
885 ! expand row 'k2' of matV directly into vecAv
886 ! non-zero range (excluding 'kn')
887 ifirst2=irangeparnz(1,k2)
888 ilast2=irangeparnz(2,k2)
889 vecav(ifirst2:ilast2)=vecav(ifirst2:ilast2)+matcoeff(ioff1+k2)* &
890 matv(ioffrow(k2)+1:ioffrow(k2)+1+ilast2-ifirst2)
891 vecav(kn2)=vecav(kn2)+matcoeff(ioff1+k2)*vecvk(k2)
892 END DO
893 ! update
894 ioff2=0
895 DO i=1,kn
896 ! correct with 2*(2v*vtAv*v^t - Av*v^t - (Av*v^t)^t)
897 DO l=max(1,i-mbnd),i
898 ioff2=ioff2+1
899 b(ioff2)=b(ioff2)+2.0_mpd*((2.0_mpd*vecn(i)*vtav-vecav(i))*vecn(l)-vecav(l)*vecn(i))
900 END DO
901 END DO
902 ! off diagonal block
903 DO i=kn+1,npar
904 ! correct with -2Av*v^t
905 DO l=max(1,i-mbnd),i
906 ioff2=ioff2+1
907 b(ioff2)=b(ioff2)-2.0_mpd*vecav(i)*vecn(l)
908 END DO
909 END DO
910 ! correct A*v for the remainung v
911 DO j2=j+1,ncon
912 k2=j2
913 IF (t) k2=ncon+1-j2
914 kn2=npar+k2-ncon
915 ioff2=int(k2-1,mpl)*int(ncon,mpl)
916 vtvp=matvtvp(ioff1+k2) ! v^t*v'
917 ! non-zero regions
918 l1=irangecoeff(1,k2)
919 l2=irangecoeff(2,k2)
920 vtavp=matvtavp(ioff2+k)
921 IF (l1 <= l2) vtavp=vtavp+dot_product(matcoeff(ioff2+l1:ioff2+l2),matvtvp(ioff1+l1:ioff1+l2)) ! v^t*(A*v')
922 l1=min(l1,k)
923 l2=max(l2,k)
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
929 END IF
930 irangecoeff(1,k2)=l1
931 irangecoeff(2,k2)=l2
932 END DO
933 END DO
934
935 CALL mpdealloc(irangecoeff)
936 CALL mpdealloc(matcoeff)
937 CALL mpdealloc(matvtavp)
938 CALL mpdealloc(matvtvp)
939 CALL mpdealloc(vecav)
940#ifdef SCOREP_USER_ENABLE
941 scorep_user_region_by_name_end("UR_qlpssq")
942#endif
943 !$POMP INST END(qlpssq)
944
945END SUBROUTINE qlpssq
946
947
955SUBROUTINE qlgete(emin,emax)
956 USE mpqldec
957
958 IMPLICIT NONE
959 INTEGER(mpi) :: i
960 INTEGER(mpi) :: ibpar
961 INTEGER(mpi) :: icoff
962 INTEGER(mpi) :: iclast
963 INTEGER(mpl) :: idiag
964
965 REAL(mpd), INTENT(OUT) :: emin
966 REAL(mpd), INTENT(OUT) :: emax
967
968 emax=matl(1)
969 emin=emax
970 DO ibpar=1,nblock ! parameter block
971 icoff=ioffblock(ibpar) ! constraint offset in parameter block
972 iclast=ioffblock(ibpar+1) ! last constraint in parameter block
973 idiag=int(ncon,mpl)*int(icoff,mpl)+1
974 DO i=icoff+1,iclast
975 IF (abs(emax) < abs(matl(idiag))) emax=matl(idiag)
976 IF (abs(emin) > abs(matl(idiag))) emin=matl(idiag)
977 idiag=idiag+ncon+1
978 END DO
979 END DO
980
981END SUBROUTINE qlgete
982
983
991SUBROUTINE qlbsub(d,y)
992 USE mpqldec
993
994 IMPLICIT NONE
995 INTEGER(mpi) :: icoff
996 INTEGER(mpi) :: iclast
997 INTEGER(mpl) :: idiag
998 INTEGER(mpi) :: k
999 INTEGER(mpi) :: nconb
1000
1001 REAL(mpd), INTENT(IN) :: d(ncon)
1002 REAL(mpd), INTENT(OUT) :: y(ncon)
1003
1004 ! solve L*y=d by forward substitution
1005 icoff=ioffblock(iblock) ! constraint offset in parameter block
1006 iclast=ioffblock(iblock+1) ! last constraint in parameter block
1007 nconb=iclast-icoff ! number of constraints in block
1008 idiag=int(ncon,mpl)*int(iclast-1,mpl)+nconb
1009 DO k=nconb,1,-1
1010 y(k)=(d(k)-dot_product(matl(idiag+1:idiag+nconb-k),y(k+1:nconb)))/matl(idiag)
1011 idiag=idiag-ncon-1
1012 END DO
1013
1014END SUBROUTINE qlbsub
1015
1018SUBROUTINE qlsetb(ib)
1019 USE mpqldec
1020
1021 IMPLICIT NONE
1022 INTEGER(mpi), INTENT(IN) :: ib
1023
1024 iblock=ib
1025
1026END SUBROUTINE qlsetb
1027
1030SUBROUTINE qldump()
1031 USE mpqldec
1032
1033 IMPLICIT NONE
1034 INTEGER(mpi) :: i
1035 INTEGER(mpi) :: ifirst
1036 INTEGER(mpi) :: ilast
1037 INTEGER(mpl) :: ioff1
1038 INTEGER(mpl) :: ioff2
1039 INTEGER(mpi) :: istat(6)
1040 INTEGER(mpi) :: j
1041 INTEGER(mpi) :: kn
1042 REAL(mpd) :: v1
1043 REAL(mpd) :: v2
1044 REAL(mpd) :: v3
1045 REAL(mpd) :: v4
1046
1047 print *
1048 ioff1=0
1049 ioff2=0
1050
1051 DO i=1, ncon
1052 kn=npar-ncon+i
1053 istat=0
1054 v1=0.;v2=0.;v3=0.;v4=0.
1055 ! expand row 'i' of matV into vecN
1056 ! non-zero range (excluding 'kn')
1057 ifirst=irangeparnz(1,i)
1058 ilast=irangeparnz(2,i)
1059 vecn=0._mpd
1060 vecn(ifirst:ilast)=matv(ioffrow(i)+1:ioffrow(i)+1+ilast-ifirst)
1061 DO j=1,npar+i-ncon
1062 IF (vecn(j) /= 0.0_mpd) THEN
1063 v2=vecn(j)
1064 IF (istat(3) == 0) THEN
1065 istat(1)=j
1066 v1=v2
1067 END IF
1068 istat(2)=j
1069 istat(3)=istat(3)+1
1070 END IF
1071 END DO
1072 ioff1=ioff1+npar
1073 DO j=1,ncon
1074 IF (matl(ioff2+j) /= 0.0_mpd) THEN
1075 v4=matl(ioff2+j)
1076 IF (istat(6) == 0) THEN
1077 istat(4)=j
1078 v3=v4
1079 END IF
1080 istat(5)=j
1081 istat(6)=istat(6)+1
1082 END IF
1083 END DO
1084 ioff2=ioff2+ncon
1085 print 100, i, istat, v1, v2, v3, v4, vecvk(i), irangeparnz(:,i)
1086 END DO
1087 print *
1088100 FORMAT(" qldump",7i8,5g13.5,2i8)
1089
1090END SUBROUTINE qldump
mpdalc::mpalloc
allocate array
Definition: mpdalc.f90:36
mpdalc::mpdealloc
deallocate array
Definition: mpdalc.f90:42
monpgs
subroutine monpgs(i)
Progress monitoring.
Definition: mpmon.f90:68
qlmrq
subroutine qlmrq(x, m, t)
Multiply right by Q(t).
Definition: mpqldec.f90:466
qldump
subroutine qldump()
Print statistics.
Definition: mpqldec.f90:1031
qlsmq
subroutine qlsmq(x, t)
Similarity transformation by Q(t).
Definition: mpqldec.f90:514
qlpssq
subroutine qlpssq(aprod, B, m, t)
Partial similarity transformation by Q(t).
Definition: mpqldec.f90:712
qldecb
subroutine qldecb(a, bpar, bcon, rcon)
QL decomposition (for disjoint block matrix).
Definition: mpqldec.f90:220
qldec
subroutine qldec(a)
QL decomposition (as single block).
Definition: mpqldec.f90:126
qlmlq
subroutine qlmlq(x, m, t)
Multiply left by Q(t) (per block).
Definition: mpqldec.f90:405
qlsetb
subroutine qlsetb(ib)
Set block.
Definition: mpqldec.f90:1019
qlbsub
subroutine qlbsub(d, y)
Backward substitution (per block).
Definition: mpqldec.f90:992
qlini
subroutine qlini(n, m, l, s, k)
Initialize QL decomposition.
Definition: mpqldec.f90:62
qlgete
subroutine qlgete(emin, emax)
Get eigenvalues.
Definition: mpqldec.f90:956
qlssq
subroutine qlssq(aprod, A, s, roff, t)
Similarity transformation by Q(t).
Definition: mpqldec.f90:574
mpdalc
(De)Allocate vectors and arrays.
Definition: mpdalc.f90:24
mpdef
Definition of constants.
Definition: mpdef.f90:24
mpdef::mpd
integer, parameter mpd
double precision
Definition: mpdef.f90:38
mpqldec
QL data.
Definition: mpqldec.f90:31
mpqldec::ncon
integer(mpi) ncon
number of constraints
Definition: mpqldec.f90:36
mpqldec::iblock
integer(mpi) iblock
active block
Definition: mpqldec.f90:39
mpqldec::vecvk
real(mpd), dimension(:), allocatable vecvk
secondary diagonal of matV (including last element)
Definition: mpqldec.f90:42
mpqldec::irangeparnz
integer(mpi), dimension(:,:), allocatable irangeparnz
range for non zero part (except vecVk)
Definition: mpqldec.f90:49
mpqldec::monpg
integer(mpi) monpg
flag for progress monitoring
Definition: mpqldec.f90:40
mpqldec::ioffrow
integer(mpl), dimension(:), allocatable ioffrow
row offsets in matV (for constrint block)
Definition: mpqldec.f90:47
mpqldec::matsize
integer(mpl) matsize
size of contraints matrix
Definition: mpqldec.f90:38
mpqldec::matv
real(mpd), dimension(:), allocatable matv
unit normals (v_i) of Householder reflectors
Definition: mpqldec.f90:41
mpqldec::matl
real(mpd), dimension(:), allocatable matl
lower diagonal matrix L
Definition: mpqldec.f90:43
mpqldec::ioffpar
integer(mpi), dimension(:), allocatable ioffpar
parameter number offsets for matV ( " )
Definition: mpqldec.f90:48
mpqldec::nblock
integer(mpi) nblock
number of blocks
Definition: mpqldec.f90:37
mpqldec::vecn
real(mpd), dimension(:), allocatable vecn
normal vector
Definition: mpqldec.f90:44
mpqldec::npar
integer(mpi) npar
number of parameters
Definition: mpqldec.f90:35
mpqldec::ioffblock
integer(mpi), dimension(:), allocatable ioffblock
block offset (1.
Definition: mpqldec.f90:46
mpqldec::nparblock
integer(mpi), dimension(:), allocatable nparblock
number of parameters in block
Definition: mpqldec.f90:45
Generated on Thu Dec 5 2024 14:18:59 for Millepede-II by doxygen 1.9.3