35 #include "qrm_common.h"
125 #if defined (sprec) || defined (dprec)
126 character,
parameter :: tran=
't', notran=
'n'
127 #elif defined (cprec) || defined (zprec)
128 character,
parameter :: tran=
'c', notran=
'n'
139 character :: uplo, unit
141 real(kind(1.d0)) :: a(lda,*)
142 real(kind(1.d0)) :: b(*)
143 integer :: mb, nb, row, col,
i, j
144 logical :: lo, ud, up
158 if(row .eq. col .and. ud)
then
161 b((j-1)*mb+
i) = (a(row,col))
165 do i= max(j-nb+n/2+1,1), mb
167 if(row .eq. col .and. ud)
then
170 b((j-1)*mb+
i) = a(row,col)
183 if(row .eq. col .and. ud)
then
186 b((j-1)*mb+
i) = a(row,col)
190 do i= max(mb-(n/2-j+1)+1,1), mb
192 if(row .eq. col .and. ud)
then
195 b((j-1)*mb+
i) = (a(row,col))
208 subroutine drpmv(uplo, trans, diag, n, a, x)
211 character :: uplo, trans, diag
213 real(kind(1.d0)) :: a(*), x(*)
216 character :: nuplo, ntrans
218 if(uplo .eq.
'l')
then
222 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
224 call dtrmv(
'u', notran, diag, n/2, a((na-n/2)*ma+1), ma, x(1), 1)
225 call dgemv(tran, ma, na-n/2-1, 1.d0, a(1), ma, &
226 & x(n/2+1), 1, 1.d0, x(1), 1)
227 call dtrmv(
'l', tran, diag, n-n/2, a((na-n/2-1)*ma+1), ma, x(n/2+1), 1)
231 call dtrmv(
'l', notran, diag, n-n/2, a((na-n/2-1)*ma+1), ma, x(n/2+1), 1)
232 call dgemv(notran, ma, na-n/2-1, 1.d0, a(1), ma, &
233 & x(1), 1, 1.d0, x(n/2+1), 1)
234 call dtrmv(
'u', tran, diag, n/2, a((na-n/2)*ma+1), ma, x(1), 1)
238 else if (uplo .eq.
'u')
then
242 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
244 call dtrmv(
'u', tran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
245 call dgemv(tran, ma-n/2-1, na, 1.d0, a(1), ma, &
246 & x(1), 1, 1.d0, x(n/2+1), 1)
247 call dtrmv(
'l', notran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
251 call dtrmv(
'l', tran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
252 call dgemv(notran, ma-n/2-1, na, 1.d0, a(1), ma, &
253 & x(n/2+1), 1, 1.d0, x(1), 1)
254 call dtrmv(
'u', notran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
268 subroutine drpsv(uplo, trans, diag, n, a, x)
271 character :: uplo, trans, diag
273 real(kind(1.d0)) :: a(*), x(*)
276 character :: nuplo, ntrans
278 if(uplo .eq.
'l')
then
282 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
284 call dtrsv(
'l', tran, diag, n-n/2, a((na-n/2-1)*ma+1), ma, x(n/2+1), 1)
285 call dgemv(tran, ma, na-n/2-1, -1.d0, a(1), ma, &
286 & x(n/2+1), 1, 1.d0, x(1), 1)
287 call dtrsv(
'u', notran, diag, n/2, a((na-n/2)*ma+1), ma, x(1), 1)
291 call dtrsv(
'u', tran, diag, n/2, a((na-n/2)*ma+1), ma, x(1), 1)
292 call dgemv(notran, ma, na-n/2-1, -1.d0, a(1), ma, &
293 & x(1), 1, 1.d0, x(n/2+1), 1)
294 call dtrsv(
'l', notran, diag, n-n/2, a((na-n/2-1)*ma+1), ma, x(n/2+1), 1)
298 else if (uplo .eq.
'u')
then
302 if((trans .eq. tran) .or. (trans .eq.
'c'))
then
304 call dtrsv(
'l', notran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
305 call dgemv(tran, ma-n/2-1, na, -1.d0, a(1), ma, &
306 & x(1), 1, 1.d0, x(n/2+1), 1)
307 call dtrsv(
'u', tran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
311 call dtrsv(
'u', notran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
312 call dgemv(notran, ma-n/2-1, na, -1.d0, a(1), ma, &
313 & x(n/2+1), 1, 1.d0, x(1), 1)
314 call dtrsv(
'l', tran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
330 subroutine drprft(direct, storev, m, k, v1, v2, ldv2, tau, t, ldt)
335 integer :: m, k, ldv2, ldt
336 character :: direct, storev
337 real(kind(1.d0)) :: v1(*), v2(ldv2,*), tau(*), t(ldt,*)
338 real(kind(1.d0)),
parameter :: zero=0.d0, one=1.d0
339 real(kind(1.d0)) :: vii, mult
340 logical,
external :: lsame
341 integer :: mv1, nv1,
i, j, off1, off2, im, in, ii, jj
349 if( lsame( direct,
'f' ) )
then
351 if( tau(
i ).eq. zero )
then
365 if( lsame( storev,
'c' ) )
then
373 #if defined (cprec) || defined (zprec)
374 if(min(im,in) .gt. 0)
then
376 call dcopy( in, v1(off2), mv1, t(ldt-in+1,1), 1 )
377 call dlacgv( in, t(ldt-in+1,1), 1)
378 call dgemv(notran, im, in, -tau(
i), &
379 & v1(off1), mv1, t(ldt-in+1,1), 1, zero, t(1,
i),1)
382 if(min(im,in) .gt. 0) call dgemv(notran, im, in, -tau(
i), &
383 & v1(off1), mv1, v1(off2), mv1, zero, t(1,
i),1)
387 off2 = max((
i-1)*mv1 + max(
i-k/2,1),1)
388 im = min(mv1,mv1-(
i-k/2)+1)
392 if(min(im,in) .gt. 0) call dgemv(tran, im, in, -tau(
i), &
393 & v1(off1), mv1, v1(off2), 1, mult, t(1,
i), 1)
395 if(min(m,
i-1) .gt. 0) call dgemv(tran, m,
i-1, -tau(
i), &
396 & v2(1,1), ldv2, v2(1,
i), 1, one, t(1,
i), 1)
399 __qrm_prnt_err(
'("_RPRFT: method not supported")')
404 call dtrmv(
'upper', notran,
'non-unit',
i-1, t, ldt, t( 1,
i ), 1 )
412 __qrm_prnt_err(
'("_RPRFT: method not supported")')
421 subroutine drprfb(side, trans, direct, storev, m, n, k, v1, v2, ldv2, t, ldt, c, ldc, work, ldwork)
425 character :: side, trans, direct, storev
426 integer :: m, n, k, ldv2, ldt, ldc, ldwork, ii
427 real(kind(1.d0)) :: v1(*), v2(ldv2,*), t(ldt,*), c(ldc,*), work(ldwork, *)
430 real(kind(1.d0)),
parameter :: one=1.d0
431 logical,
external :: lsame
434 if( m.le.0 .or. n.le.0 )
return
436 if( lsame( trans,
'n' ) )
then
437 #if defined (sprec) || defined (dprec)
439 #elif defined (cprec) || defined (zprec)
446 if( lsame( storev,
'c' ) )
then
447 if( lsame( direct,
'f' ) )
then
451 if( lsame( side,
'l' ) )
then
458 call dcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
459 #if defined(cprec) || defined(zprec)
460 call dlacgv(n, work(1,j), 1)
464 call
drpmm(
'right',
'lower', notran,
'unit', n, k, one, v1, work, ldwork )
467 call dgemm( tran, notran, n, k, m-k, one, c( k+1, 1 ), ldc, &
468 & v2( 1, 1 ), ldv2, one, work, ldwork )
472 call dtrmm(
'right',
'upper', transt,
'non-unit', n, k, one, t, ldt, work, ldwork )
476 call dgemm( notran, tran, m-k, n, k, -one, v2( 1, 1 ), ldv2, work, &
477 & ldwork, one, c( k+1, 1 ), ldc )
480 call
drpmm(
'right',
'lower', tran,
'unit', n, k, one, v1, work, ldwork )
485 c( j,
i ) = c( j,
i ) - (work(
i, j ))
489 else if( lsame( side,
'r' ) )
then
490 __qrm_prnt_err(
'("_RPRFB: method not supported")')
493 __qrm_prnt_err(
'("_RPRFB: method not supported")')
497 __qrm_prnt_err(
'("_RPRFB: method not supported")')
507 subroutine drpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
510 character :: uplo, trans, diag, side
512 real(kind(1.d0)) :: a(*), x(ldx,*), alpha
514 integer :: ma, na, p1, p2, p3, off
515 logical,
external :: lsame
517 if (lsame(side,
'l'))
then
518 if(uplo .eq.
'l')
then
526 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
528 call dtrmm(
'l',
'u', notran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
529 call dgemm(tran, notran, m/2, n, ma, alpha, a(p1), ma, x(off,1), ldx, 1.d0, x(1,1), ldx)
530 call dtrmm(
'l',
'l', tran, diag, ma, n, alpha, a(p2), ma, x(off,1), ldx)
534 call dtrmm(
'l',
'l', notran, diag, ma, n, alpha, a(p2), ma, x(off,1), ldx)
535 call dgemm(notran, notran, ma, n, m/2, alpha, a(p1), ma, x(1,1), ldx, 1.d0, x(off,1), ldx)
536 call dtrmm(
'l',
'u', tran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
540 else if (uplo .eq.
'u')
then
548 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
550 call dtrmm(
'l',
'u', tran, diag, na, n, alpha, a(p2), ma, x(off,1), ldx)
551 call dgemm(tran, notran, na, n, m/2, alpha, a(p1), ma, x(1,1), ldx, 1.d0, x(off,1), ldx)
552 call dtrmm(
'l',
'l', notran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
556 call dtrmm(
'l',
'l', tran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
557 call dgemm(notran, notran, m/2, n, na, alpha, a(p1), ma, x(off,1), ldx, 1.d0, x(1,1), ldx)
558 call dtrmm(
'l',
'u', notran, diag, na, n, alpha, a(p2), ma, x(off, 1), ldx)
566 if(uplo .eq.
'l')
then
574 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
576 call dtrmm(
'r',
'l', tran, diag, m, ma, alpha, a(p2), ma, x(1,off), ldx)
577 call dgemm(notran, tran, m, ma, n/2, alpha, x(1,1), ldx, a(p1), ma, 1.d0, x(1,off), ldx)
578 call dtrmm(
'r',
'u', notran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
581 call dtrmm(
'r',
'u', tran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
582 call dgemm(notran, notran, m, n/2, ma, alpha, x(1,off), ldx, a(p1), ma, 1.d0, x(1,1), ldx)
583 call dtrmm(
'r',
'l', notran, diag, m, ma, alpha, a(p2), ma, x(1,off), ldx)
587 else if (uplo .eq.
'u')
then
595 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
597 call dtrmm(
'r',
'l', notran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
598 call dgemm(notran, tran, m, n/2, na, alpha, x(1,off), ldx, a(p1), ma, 1.d0, x(1,1), ldx)
599 call dtrmm(
'r',
'u', tran, diag, m, na, alpha, a(p2), ma, x(1,off), ldx)
603 call dtrmm(
'r',
'u', notran, diag, m, na, alpha, a(p2), ma, x(1, off), ldx)
604 call dgemm(notran, notran, m, na, n/2, alpha, x(1,1), ldx, a(p1), ma, 1.d0, x(1,off), ldx)
605 call dtrmm(
'r',
'l', tran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
621 subroutine drpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
624 character :: uplo, trans, diag, side
626 real(kind(1.d0)) :: a(*), x(ldx,*), alpha
628 integer :: ma, na, p1, p2, p3, off, ii
629 logical,
external :: lsame
631 if (lsame(side,
'l'))
then
632 if(uplo .eq.
'l')
then
640 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
642 call dtrsm(
'l',
'l', tran, diag, ma, n, alpha, a(p2), ma, x(off,1), ldx)
643 call dgemm(tran, notran, m/2, n, ma, -1.d0, a(p1), ma, x(off,1), ldx, alpha, x(1,1), ldx)
644 call dtrsm(
'l',
'u', notran, diag, m/2, n, 1.d0, a(p3), ma, x(1,1), ldx)
648 call dtrsm(
'l',
'u', tran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
649 call dgemm(notran, notran, ma, n, m/2, -1.d0, a(p1), ma, x(1,1), ldx, alpha, x(off,1), ldx)
650 call dtrsm(
'l',
'l', notran, diag, ma, n, 1.d0, a(p2), ma, x(off,1), ldx)
654 else if (uplo .eq.
'u')
then
662 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
664 call dtrsm(
'l',
'l', notran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
665 call dgemm(tran, notran, na, n, m/2, -1.d0, a(p1), ma, x(1,1), ldx, alpha, x(off,1), ldx)
666 call dtrsm(
'l',
'u', tran, diag, na, n, 1.d0, a(p2), ma, x(off,1), ldx)
670 call dtrsm(
'l',
'u', notran, diag, na, n, alpha, a(p2), ma, x(off, 1), ldx)
671 call dgemm(notran, notran, m/2, n, na, -1.d0, a(p1), ma, x(off,1), ldx, alpha, x(1,1), ldx)
672 call dtrsm(
'l',
'l', tran, diag, m/2, n, 1.d0, a(p3), ma, x(1,1), ldx)
680 if(uplo .eq.
'l')
then
688 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
690 call dtrsm(
'r',
'u', notran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
691 call dgemm(notran, tran, m, ma, n/2, -1.d0, x(1,1), ldx, a(p1), ma, alpha, x(1,off), ldx)
692 call dtrsm(
'r',
'l', tran, diag, m, ma, 1.d0, a(p2), ma, x(1,off), ldx)
696 call dtrsm(
'r',
'l', notran, diag, m, ma, alpha, a(p2), ma, x(1,off), ldx)
697 call dgemm(notran, notran, m, n/2, ma, -1.d0, x(1,off), ldx, a(p1), ma, alpha, x(1,1), ldx)
698 call dtrsm(
'r',
'u', tran, diag, m, n/2, 1.d0, a(p3), ma, x(1,1), ldx)
702 else if (uplo .eq.
'u')
then
710 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
712 call dtrsm(
'r',
'u', tran, diag, m, na, alpha, a(p2), ma, x(1,off), ldx)
713 call dgemm(notran, tran, m, n/2, na, -1.d0, x(1,off), ldx, a(p1), ma, alpha, x(1,1), ldx)
714 call dtrsm(
'r',
'l', notran, diag, m, n/2, 1.d0, a(p3), ma, x(1,1), ldx)
718 call dtrsm(
'r',
'l', tran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
719 call dgemm(notran, notran, m, na, n/2, -1.d0, x(1,1), ldx, a(p1), ma, alpha, x(1,off), ldx)
720 call dtrsm(
'r',
'u', notran, diag, m, na, 1.d0, a(p2), ma, x(1, off), ldx)
737 real(kind(1.d0)) :: b(*)
739 integer ::
i, j, mb, nb
741 if(uplo .eq.
'l')
then
744 else if (uplo .eq.
'u')
then
752 write(*,
'(f6.1, x)', advance=
'no')b((j-1)*mb +
i)
This module contains all the error management routines and data.
subroutine dqrm_print_rfpf(uplo, n, b)
subroutine drpsv(uplo, trans, diag, n, a, x)
subroutine drpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
subroutine drprfb(side, trans, direct, storev, m, n, k, v1, v2, ldv2, t, ldt, c, ldc, work, ldwork)
subroutine drprft(direct, storev, m, k, v1, v2, ldv2, tau, t, ldt)
subroutine drpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
subroutine dqrm_to_rfpf(uplo, unit, n, a, lda, b)
subroutine drpmv(uplo, trans, diag, n, a, x)
This module contains various string handling routines.
This module contains an implementation of some operations on triangular/trapezoidal matrices stored i...