QR_MUMPS
 All Classes Files Functions Variables Enumerations Enumerator Pages
qrm_rfpf_mod.F90
Go to the documentation of this file.
1 !! ##############################################################################################
2 !!
3 !! Copyright 2012 CNRS, INPT
4 !!
5 !! This file is part of qr_mumps.
6 !!
7 !! qr_mumps is free software: you can redistribute it and/or modify
8 !! it under the terms of the GNU Lesser General Public License as
9 !! published by the Free Software Foundation, either version 3 of
10 !! the License, or (at your option) any later version.
11 !!
12 !! qr_mumps is distributed in the hope that it will be useful,
13 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
14 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 !! GNU Lesser General Public License for more details.
16 !!
17 !! You can find a copy of the GNU Lesser General Public License
18 !! in the qr_mumps/doc directory.
19 !!
20 !! ##############################################################################################
21 
22 
23 !! ##############################################################################################
33 
34 
35 #include "qrm_common.h"
36 
39 
122 
123  use qrm_error_mod
124 
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'
129 #endif
130 
131 
132 contains
133 
134  subroutine _qrm_to_rfpf(uplo, unit, n, a, lda, b)
135 
136  use qrm_string_mod
137  implicit none
138 
139  character :: uplo, unit
140  integer :: n, lda
141  _qrm_data :: a(lda,*)
142  _qrm_data :: b(*)
143  integer :: mb, nb, row, col, i, j
144  logical :: lo, ud, up
145 
146  lo = qrm_str_tolower(uplo) .eq. 'l'
147  up = qrm_str_tolower(uplo) .eq. 'u'
148  ud = qrm_str_tolower(unit) .eq. 'u'
149 
150  if(lo) then
151  mb = (n+1)/2
152  nb = n+(n/2-mb+1)
153 
154  do j=1, nb
155  row = j-nb+n/2
156  do i=1, j-nb+n/2
157  col = i
158  if(row .eq. col .and. ud) then
159  b((j-1)*mb+i) = _qrm_one
160  else
161  b((j-1)*mb+i) = _conjg(a(row,col))
162  end if
163  end do
164  col = j
165  do i= max(j-nb+n/2+1,1), mb
166  row = i+n/2
167  if(row .eq. col .and. ud) then
168  b((j-1)*mb+i) = _qrm_one
169  else
170  b((j-1)*mb+i) = a(row,col)
171  end if
172  end do
173  end do
174 
175  else if (up) then
176  nb = (n+1)/2
177  mb = n+(n/2-nb+1)
178 
179  do j=1, nb
180  col = j+n/2
181  do i=1, mb-(n/2-j+1)
182  row = i
183  if(row .eq. col .and. ud) then
184  b((j-1)*mb+i) = _qrm_one
185  else
186  b((j-1)*mb+i) = a(row,col)
187  end if
188  end do
189  row = j
190  do i= max(mb-(n/2-j+1)+1,1), mb
191  col = i - (mb-n/2)
192  if(row .eq. col .and. ud) then
193  b((j-1)*mb+i) = _qrm_one
194  else
195  b((j-1)*mb+i) = _conjg(a(row,col))
196  end if
197  end do
198  end do
199 
200  end if
201 
202  return
203 
204  end subroutine _qrm_to_rfpf
205 
206 
207 
208  subroutine _xrpmv(uplo, trans, diag, n, a, x)
209 
210  implicit none
211  character :: uplo, trans, diag
212  integer :: n
213  _qrm_data :: a(*), x(*)
214 
215  integer :: ma, na
216  character :: nuplo, ntrans
217 
218  if(uplo .eq. 'l') then
219  ma = (n+1)/2
220  na = n+(n/2-ma+1)
221 
222  if((trans .eq. 't') .or. (trans .eq. 'c')) then
223 
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)
228 
229  else
230 
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)
235 
236  end if
237 
238  else if (uplo .eq. 'u') then
239  na = (n+1)/2
240  ma = n+(n/2-na+1)
241 
242  if((trans .eq. 't') .or. (trans .eq. 'c')) then
243 
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)
248 
249  else
250 
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)
255 
256  end if
257 
258  end if
259 
260  return
261 
262  end subroutine _xrpmv
263 
264 
265 
266 
267 
268  subroutine _xrpsv(uplo, trans, diag, n, a, x)
269 
270  implicit none
271  character :: uplo, trans, diag
272  integer :: n
273  _qrm_data :: a(*), x(*)
274 
275  integer :: ma, na
276  character :: nuplo, ntrans
277 
278  if(uplo .eq. 'l') then
279  ma = (n+1)/2
280  na = n+(n/2-ma+1)
281 
282  if((trans .eq. 't') .or. (trans .eq. 'c')) then
283 
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)
288 
289  else
290 
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)
295 
296  end if
297 
298  else if (uplo .eq. 'u') then
299  na = (n+1)/2
300  ma = n+(n/2-na+1)
301 
302  if((trans .eq. tran) .or. (trans .eq. 'c')) then
303 
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)
308 
309  else
310 
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)
315 
316  end if
317 
318  end if
319 
320  return
321 
322  end subroutine _xrpsv
323 
324 
325 
326 
327 
328 
329 
330  subroutine _xrprft(direct, storev, m, k, v1, v2, ldv2, tau, t, ldt)
331  ! This routine is different from LAPACK's one because the V matrix
332  ! already has ones on the diagonal.
333  use qrm_const_mod
334  implicit none
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
342  character :: tr
343 
344  if( k.eq.0 ) return
345 
346  mv1 = (k+1)/2
347  nv1 = k+(k/2-mv1+1)
348 
349  if( lsame( direct, 'f' ) ) then
350  do i = 1, k
351  if( tau( i ).eq. zero ) then
352  !
353  ! h(i) = i
354  !
355  do j = 1, i
356  t( j, i ) = zero
357  end do
358  else
359  if(i .le. k/2) then
360  mult = _qrm_one
361  else
362  mult = _qrm_zero
363  end if
364 
365  if( lsame( storev, 'c' ) ) then
366  off1 = (k/2+i)*mv1+1
367  off2 = (k/2+i)*mv1+i
368  im = i-1
369  in = k/2-i+1
370 
371 
372  ! FIXME: check this (seems to work). last good revision 293 (non tread safe)
373 #if defined (cprec) || defined (zprec)
374  if(min(im,in) .gt. 0) then
375  ! using the lower part ot T as a scratch memory
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)
380  end if
381 #else
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)
384 #endif
385 
386  off1 = max(i-k/2,1)
387  off2 = max((i-1)*mv1 + max(i-k/2,1),1)
388  im = min(mv1,mv1-(i-k/2)+1)
389  in = i-1
390 
391 
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)
394 
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)
397 
398  else
399  __qrm_prnt_err('("_RPRFT: method not supported")')
400  return
401  end if
402 
403 
404  call _xtrmv( 'upper', notran, 'non-unit', i-1, t, ldt, t( 1, i ), 1 )
405 
406  t( i, i ) = tau( i )
407 
408 
409  end if
410  end do
411  else
412  __qrm_prnt_err('("_RPRFT: method not supported")')
413  return
414  end if
415  return
416 
417  end subroutine _xrprft
418 
419 
420 
421  subroutine _xrprfb(side, trans, direct, storev, m, n, k, v1, v2, ldv2, t, ldt, c, ldc, work, ldwork)
422  use qrm_const_mod
423  implicit none
424 
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, *)
428 
429  integer :: i, j
430  _qrm_data, parameter :: one=_qrm_one
431  logical, external :: lsame
432  character :: transt
433 
434  if( m.le.0 .or. n.le.0 ) return
435 
436  if( lsame( trans, 'n' ) ) then
437 #if defined (sprec) || defined (dprec)
438  transt = 't'
439 #elif defined (cprec) || defined (zprec)
440  transt = 'c'
441 #endif
442  else
443  transt = 'n'
444  end if
445 
446  if( lsame( storev, 'c' ) ) then
447  if( lsame( direct, 'f' ) ) then
448  ! let v = ( v1 ) (first k rows)
449  ! ( v2 )
450  ! where v1 is unit lower triangular.
451  if( lsame( side, 'l' ) ) then
452  ! form h * c or h' * c where c = ( c1 )
453  ! ( c2 )
454  !
455  ! w := c' * v = (c1'*v1 + c2'*v2) (stored in work)
456  ! w := c1'
457  do j = 1, k
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)
461 #endif
462  end do
463  ! w := w * v1
464  call _xrpmm( 'right', 'lower', notran, 'unit', n, k, one, v1, work, ldwork )
465  if( m.gt.k ) then
466  ! w := w + c2'*v2
467  call _xgemm( tran, notran, n, k, m-k, one, c( k+1, 1 ), ldc, &
468  & v2( 1, 1 ), ldv2, one, work, ldwork )
469  end if
470  ! w := w * t' or w * t
471  !FIXME not thread-safe (for some reason)
472  call _xtrmm( 'right', 'upper', transt, 'non-unit', n, k, one, t, ldt, work, ldwork )
473  ! c := c - v * w'
474  if( m.gt.k ) then
475  ! c2 := c2 - v2 * w'
476  call _xgemm( notran, tran, m-k, n, k, -one, v2( 1, 1 ), ldv2, work, &
477  & ldwork, one, c( k+1, 1 ), ldc )
478  end if
479  ! w := w * v1'
480  call _xrpmm( 'right', 'lower', tran, 'unit', n, k, one, v1, work, ldwork )
481 
482  ! c1 := c1 - w'
483  do j = 1, k
484  do i = 1, n
485  c( j, i ) = c( j, i ) - _conjg(work( i, j ))
486  end do
487  end do
488 
489  else if( lsame( side, 'r' ) ) then
490  __qrm_prnt_err('("_RPRFB: method not supported")')
491  end if
492  else
493  __qrm_prnt_err('("_RPRFB: method not supported")')
494  end if
495 
496  else
497  __qrm_prnt_err('("_RPRFB: method not supported")')
498  end if
499 
500  return
501  end subroutine _xrprfb
502 
503 
504 
505 
506 
507  subroutine _xrpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
508 
509  implicit none
510  character :: uplo, trans, diag, side
511  integer :: m, n, ldx
512  _qrm_data :: a(*), x(ldx,*), alpha
513 
514  integer :: ma, na, p1, p2, p3, off
515  logical, external :: lsame
516 
517  if (lsame(side, 'l')) then
518  if(uplo .eq. 'l') then
519  ma = (m+1)/2
520  na = m+(m/2-ma+1)
521  p1 = 1
522  p2 = (m/2)*ma+1
523  p3 = (m/2+1)*ma+1
524  off = m/2+1
525 
526  if((trans .eq. 't') .or. (trans .eq. 'c')) then
527 
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)
531 
532  else
533 
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)
537 
538  end if
539 
540  else if (uplo .eq. 'u') then
541  na = (m+1)/2
542  ma = m+(m/2-na+1)
543  p1 = 1
544  p2 = m/2+1
545  p3 = m/2+2
546  off = m/2+1
547 
548  if((trans .eq. 't') .or. (trans .eq. 'c')) then
549 
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)
553 
554  else
555 
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)
559 
560  end if
561 
562  end if
563 
564  else
565 
566  if(uplo .eq. 'l') then
567  ma = (n+1)/2
568  na = n+(n/2-ma+1)
569  p1 = 1
570  p2 = (n/2)*ma+1
571  p3 = (n/2+1)*ma+1
572  off = n/2+1
573 
574  if((trans .eq. 't') .or. (trans .eq. 'c')) then
575 
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)
579 
580  else
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)
584 
585  end if
586 
587  else if (uplo .eq. 'u') then
588  na = (n+1)/2
589  ma = n+(n/2-na+1)
590  p1 = 1
591  p2 = n/2+1
592  p3 = n/2+2
593  off = n/2+1
594 
595  if((trans .eq. 't') .or. (trans .eq. 'c')) then
596 
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)
600 
601  else
602 
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)
606 
607  end if
608 
609  end if
610 
611  end if
612 
613 
614 
615  return
616 
617  end subroutine _xrpmm
618 
619 
620 
621  subroutine _xrpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
622 
623  implicit none
624  character :: uplo, trans, diag, side
625  integer :: m, n, ldx
626  _qrm_data :: a(*), x(ldx,*), alpha
627 
628  integer :: ma, na, p1, p2, p3, off, ii
629  logical, external :: lsame
630 
631  if (lsame(side, 'l')) then
632  if(uplo .eq. 'l') then
633  ma = (m+1)/2
634  na = m+(m/2-ma+1)
635  p1 = 1
636  p2 = (m/2)*ma+1
637  p3 = (m/2+1)*ma+1
638  off = m/2+1
639 
640  if((trans .eq. 't') .or. (trans .eq. 'c')) then
641 
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)
645 
646  else
647 
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)
651 
652  end if
653 
654  else if (uplo .eq. 'u') then
655  na = (m+1)/2
656  ma = m+(m/2-na+1)
657  p1 = 1
658  p2 = m/2+1
659  p3 = m/2+2
660  off = m/2+1
661 
662  if((trans .eq. 't') .or. (trans .eq. 'c')) then
663 
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)
667 
668  else
669 
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)
673 
674  end if
675 
676  end if
677 
678  else
679 
680  if(uplo .eq. 'l') then
681  ma = (n+1)/2
682  na = n+(n/2-ma+1)
683  p1 = 1
684  p2 = (n/2)*ma+1
685  p3 = (n/2+1)*ma+1
686  off = n/2+1
687 
688  if((trans .eq. 't') .or. (trans .eq. 'c')) then
689 
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)
693 
694  else
695 
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)
699 
700  end if
701 
702  else if (uplo .eq. 'u') then
703  na = (n+1)/2
704  ma = n+(n/2-na+1)
705  p1 = 1
706  p2 = n/2+1
707  p3 = n/2+2
708  off = n/2+1
709 
710  if((trans .eq. 't') .or. (trans .eq. 'c')) then
711 
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)
715 
716  else
717 
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)
721 
722  end if
723 
724  end if
725 
726  end if
727 
728  return
729 
730  end subroutine _xrpsm
731 
732 
733  subroutine _qrm_print_rfpf(uplo, n, b)
734 
735  character :: uplo
736  integer :: n
737  _qrm_data :: b(*)
738 
739  integer :: i, j, mb, nb
740 
741  if(uplo .eq. 'l') then
742  mb = (n+1)/2
743  nb = n+(n/2-mb+1)
744  else if (uplo .eq. 'u') then
745  nb = (n+1)/2
746  mb = n+(n/2-nb+1)
747  end if
748 
749 
750  do i=1, mb
751  do j=1, nb
752  write(*,'(f6.1, x)', advance='no')b((j-1)*mb + i)
753  end do
754  write(*,'(" ")')
755  end do
756 
757  end subroutine _qrm_print_rfpf
758 
759 
760 
761 end module _qrm_rfpf_mod
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.
int i
Definition: secs.c:40
subroutine _xrpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
subroutine _xrpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)