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 _qrm_data :: a(lda,*)
143 integer :: mb, nb, row, col,
i, j
144 logical :: lo, ud, up
158 if(row .eq. col .and. ud)
then
159 b((j-1)*mb+
i) = _qrm_one
161 b((j-1)*mb+
i) = _conjg(a(row,col))
165 do i= max(j-nb+n/2+1,1), mb
167 if(row .eq. col .and. ud)
then
168 b((j-1)*mb+
i) = _qrm_one
170 b((j-1)*mb+
i) = a(row,col)
183 if(row .eq. col .and. ud)
then
184 b((j-1)*mb+
i) = _qrm_one
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
193 b((j-1)*mb+
i) = _qrm_one
195 b((j-1)*mb+
i) = _conjg(a(row,col))
208 subroutine _xrpmv(uplo, trans, diag, n, a, x)
211 character :: uplo, trans, diag
213 _qrm_data :: a(*), x(*)
216 character :: nuplo, ntrans
218 if(uplo .eq.
'l')
then
222 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
224 call _xtrmv(
'u', notran, diag, n/2, a((na-n/2)*ma+1), ma, x(1), 1)
225 call _xgemv(tran, ma, na-n/2-1, _qrm_one, a(1), ma, &
226 & x(n/2+1), 1, _qrm_one, x(1), 1)
227 call _xtrmv(
'l', tran, diag, n-n/2, a((na-n/2-1)*ma+1), ma, x(n/2+1), 1)
231 call _xtrmv(
'l', notran, diag, n-n/2, a((na-n/2-1)*ma+1), ma, x(n/2+1), 1)
232 call _xgemv(notran, ma, na-n/2-1, _qrm_one, a(1), ma, &
233 & x(1), 1, _qrm_one, x(n/2+1), 1)
234 call _xtrmv(
'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 _xtrmv(
'u', tran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
245 call _xgemv(tran, ma-n/2-1, na, _qrm_one, a(1), ma, &
246 & x(1), 1, _qrm_one, x(n/2+1), 1)
247 call _xtrmv(
'l', notran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
251 call _xtrmv(
'l', tran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
252 call _xgemv(notran, ma-n/2-1, na, _qrm_one, a(1), ma, &
253 & x(n/2+1), 1, _qrm_one, x(1), 1)
254 call _xtrmv(
'u', notran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
268 subroutine _xrpsv(uplo, trans, diag, n, a, x)
271 character :: uplo, trans, diag
273 _qrm_data :: a(*), x(*)
276 character :: nuplo, ntrans
278 if(uplo .eq.
'l')
then
282 if((trans .eq.
't') .or. (trans .eq.
'c'))
then
284 call _xtrsv(
'l', tran, diag, n-n/2, a((na-n/2-1)*ma+1), ma, x(n/2+1), 1)
285 call _xgemv(tran, ma, na-n/2-1, -_qrm_one, a(1), ma, &
286 & x(n/2+1), 1, _qrm_one, x(1), 1)
287 call _xtrsv(
'u', notran, diag, n/2, a((na-n/2)*ma+1), ma, x(1), 1)
291 call _xtrsv(
'u', tran, diag, n/2, a((na-n/2)*ma+1), ma, x(1), 1)
292 call _xgemv(notran, ma, na-n/2-1, -_qrm_one, a(1), ma, &
293 & x(1), 1, _qrm_one, x(n/2+1), 1)
294 call _xtrsv(
'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 _xtrsv(
'l', notran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
305 call _xgemv(tran, ma-n/2-1, na, -_qrm_one, a(1), ma, &
306 & x(1), 1, _qrm_one, x(n/2+1), 1)
307 call _xtrsv(
'u', tran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
311 call _xtrsv(
'u', notran, diag, na, a((ma-n/2-1)+1), ma, x(n/2+1), 1)
312 call _xgemv(notran, ma-n/2-1, na, -_qrm_one, a(1), ma, &
313 & x(n/2+1), 1, _qrm_one, x(1), 1)
314 call _xtrsv(
'l', tran, diag, n/2, a((ma-n/2)+1), ma, x(1), 1)
330 subroutine _xrprft(direct, storev, m, k, v1, v2, ldv2, tau, t, ldt)
335 integer :: m, k, ldv2, ldt
336 character :: direct, storev
337 _qrm_data :: v1(*), v2(ldv2,*), tau(*), t(ldt,*)
338 _qrm_data,
parameter :: zero=_qrm_zero, one=_qrm_one
339 _qrm_data :: 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 _xcopy( in, v1(off2), mv1, t(ldt-in+1,1), 1 )
377 call _xlacgv( in, t(ldt-in+1,1), 1)
378 call _xgemv(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 _xgemv(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 _xgemv(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 _xgemv(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 _xtrmv(
'upper', notran,
'non-unit',
i-1, t, ldt, t( 1,
i ), 1 )
412 __qrm_prnt_err(
'("_RPRFT: method not supported")')
421 subroutine _xrprfb(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 _qrm_data :: v1(*), v2(ldv2,*), t(ldt,*), c(ldc,*), work(ldwork, *)
430 _qrm_data,
parameter :: one=_qrm_one
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 _xcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
459 #if defined(cprec) || defined(zprec)
460 call _xlacgv(n, work(1,j), 1)
464 call
_xrpmm(
'right',
'lower', notran,
'unit', n, k, one, v1, work, ldwork )
467 call _xgemm( tran, notran, n, k, m-k, one, c( k+1, 1 ), ldc, &
468 & v2( 1, 1 ), ldv2, one, work, ldwork )
472 call _xtrmm(
'right',
'upper', transt,
'non-unit', n, k, one, t, ldt, work, ldwork )
476 call _xgemm( notran, tran, m-k, n, k, -one, v2( 1, 1 ), ldv2, work, &
477 & ldwork, one, c( k+1, 1 ), ldc )
480 call
_xrpmm(
'right',
'lower', tran,
'unit', n, k, one, v1, work, ldwork )
485 c( j,
i ) = c( j,
i ) - _conjg(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 _xrpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
510 character :: uplo, trans, diag, side
512 _qrm_data :: 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 _xtrmm(
'l',
'u', notran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
529 call _xgemm(tran, notran, m/2, n, ma, alpha, a(p1), ma, x(off,1), ldx, _qrm_one, x(1,1), ldx)
530 call _xtrmm(
'l',
'l', tran, diag, ma, n, alpha, a(p2), ma, x(off,1), ldx)
534 call _xtrmm(
'l',
'l', notran, diag, ma, n, alpha, a(p2), ma, x(off,1), ldx)
535 call _xgemm(notran, notran, ma, n, m/2, alpha, a(p1), ma, x(1,1), ldx, _qrm_one, x(off,1), ldx)
536 call _xtrmm(
'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 _xtrmm(
'l',
'u', tran, diag, na, n, alpha, a(p2), ma, x(off,1), ldx)
551 call _xgemm(tran, notran, na, n, m/2, alpha, a(p1), ma, x(1,1), ldx, _qrm_one, x(off,1), ldx)
552 call _xtrmm(
'l',
'l', notran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
556 call _xtrmm(
'l',
'l', tran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
557 call _xgemm(notran, notran, m/2, n, na, alpha, a(p1), ma, x(off,1), ldx, _qrm_one, x(1,1), ldx)
558 call _xtrmm(
'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 _xtrmm(
'r',
'l', tran, diag, m, ma, alpha, a(p2), ma, x(1,off), ldx)
577 call _xgemm(notran, tran, m, ma, n/2, alpha, x(1,1), ldx, a(p1), ma, _qrm_one, x(1,off), ldx)
578 call _xtrmm(
'r',
'u', notran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
581 call _xtrmm(
'r',
'u', tran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
582 call _xgemm(notran, notran, m, n/2, ma, alpha, x(1,off), ldx, a(p1), ma, _qrm_one, x(1,1), ldx)
583 call _xtrmm(
'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 _xtrmm(
'r',
'l', notran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
598 call _xgemm(notran, tran, m, n/2, na, alpha, x(1,off), ldx, a(p1), ma, _qrm_one, x(1,1), ldx)
599 call _xtrmm(
'r',
'u', tran, diag, m, na, alpha, a(p2), ma, x(1,off), ldx)
603 call _xtrmm(
'r',
'u', notran, diag, m, na, alpha, a(p2), ma, x(1, off), ldx)
604 call _xgemm(notran, notran, m, na, n/2, alpha, x(1,1), ldx, a(p1), ma, _qrm_one, x(1,off), ldx)
605 call _xtrmm(
'r',
'l', tran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
621 subroutine _xrpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
624 character :: uplo, trans, diag, side
626 _qrm_data :: 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 _xtrsm(
'l',
'l', tran, diag, ma, n, alpha, a(p2), ma, x(off,1), ldx)
643 call _xgemm(tran, notran, m/2, n, ma, -_qrm_one, a(p1), ma, x(off,1), ldx, alpha, x(1,1), ldx)
644 call _xtrsm(
'l',
'u', notran, diag, m/2, n, _qrm_one, a(p3), ma, x(1,1), ldx)
648 call _xtrsm(
'l',
'u', tran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
649 call _xgemm(notran, notran, ma, n, m/2, -_qrm_one, a(p1), ma, x(1,1), ldx, alpha, x(off,1), ldx)
650 call _xtrsm(
'l',
'l', notran, diag, ma, n, _qrm_one, 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 _xtrsm(
'l',
'l', notran, diag, m/2, n, alpha, a(p3), ma, x(1,1), ldx)
665 call _xgemm(tran, notran, na, n, m/2, -_qrm_one, a(p1), ma, x(1,1), ldx, alpha, x(off,1), ldx)
666 call _xtrsm(
'l',
'u', tran, diag, na, n, _qrm_one, a(p2), ma, x(off,1), ldx)
670 call _xtrsm(
'l',
'u', notran, diag, na, n, alpha, a(p2), ma, x(off, 1), ldx)
671 call _xgemm(notran, notran, m/2, n, na, -_qrm_one, a(p1), ma, x(off,1), ldx, alpha, x(1,1), ldx)
672 call _xtrsm(
'l',
'l', tran, diag, m/2, n, _qrm_one, 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 _xtrsm(
'r',
'u', notran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
691 call _xgemm(notran, tran, m, ma, n/2, -_qrm_one, x(1,1), ldx, a(p1), ma, alpha, x(1,off), ldx)
692 call _xtrsm(
'r',
'l', tran, diag, m, ma, _qrm_one, a(p2), ma, x(1,off), ldx)
696 call _xtrsm(
'r',
'l', notran, diag, m, ma, alpha, a(p2), ma, x(1,off), ldx)
697 call _xgemm(notran, notran, m, n/2, ma, -_qrm_one, x(1,off), ldx, a(p1), ma, alpha, x(1,1), ldx)
698 call _xtrsm(
'r',
'u', tran, diag, m, n/2, _qrm_one, 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 _xtrsm(
'r',
'u', tran, diag, m, na, alpha, a(p2), ma, x(1,off), ldx)
713 call _xgemm(notran, tran, m, n/2, na, -_qrm_one, x(1,off), ldx, a(p1), ma, alpha, x(1,1), ldx)
714 call _xtrsm(
'r',
'l', notran, diag, m, n/2, _qrm_one, a(p3), ma, x(1,1), ldx)
718 call _xtrsm(
'r',
'l', tran, diag, m, n/2, alpha, a(p3), ma, x(1,1), ldx)
719 call _xgemm(notran, notran, m, na, n/2, -_qrm_one, x(1,1), ldx, a(p1), ma, alpha, x(1,off), ldx)
720 call _xtrsm(
'r',
'u', notran, diag, m, na, _qrm_one, a(p2), ma, x(1, off), ldx)
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.
This module contains an implementation of some operations on triangular/trapezoidal matrices stored i...
subroutine _xrprfb(side, trans, direct, storev, m, n, k, v1, v2, ldv2, t, ldt, c, ldc, work, ldwork)
subroutine _qrm_print_rfpf(uplo, n, b)
subroutine _qrm_to_rfpf(uplo, unit, n, a, lda, b)
subroutine _xrpmv(uplo, trans, diag, n, a, x)
subroutine _xrpsv(uplo, trans, diag, n, a, x)
subroutine _xrprft(direct, storev, m, k, v1, v2, ldv2, tau, t, ldt)
This module contains various string handling routines.
subroutine _xrpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
subroutine _xrpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)