QR_MUMPS
dqrm_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 dqrm_to_rfpf(uplo, unit, n, a, lda, b)
136  use qrm_string_mod
137  implicit none
138 
139  character :: uplo, unit
140  integer :: n, lda
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
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) = 1.d0
160  else
161  b((j-1)*mb+i) = (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) = 1.d0
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) = 1.d0
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) = 1.d0
194  else
195  b((j-1)*mb+i) = (a(row,col))
196  end if
197  end do
198  end do
199 
200  end if
201 
202  return
203 
204  end subroutine dqrm_to_rfpf
205 
206 
207 
208  subroutine drpmv(uplo, trans, diag, n, a, x)
210  implicit none
211  character :: uplo, trans, diag
212  integer :: n
213  real(kind(1.d0)) :: 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 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)
228 
229  else
230 
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)
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 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)
248 
249  else
250 
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)
255 
256  end if
257 
258  end if
259 
260  return
261 
262  end subroutine drpmv
263 
264 
265 
266 
267 
268  subroutine drpsv(uplo, trans, diag, n, a, x)
270  implicit none
271  character :: uplo, trans, diag
272  integer :: n
273  real(kind(1.d0)) :: 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 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)
288 
289  else
290 
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)
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 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)
308 
309  else
310 
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)
315 
316  end if
317 
318  end if
319 
320  return
321 
322  end subroutine drpsv
323 
324 
325 
326 
327 
328 
329 
330  subroutine drprft(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  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
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 = 1.d0
361  else
362  mult = 0.d0
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 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)
380  end if
381 #else
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)
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 dgemv(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 dgemv(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 dtrmv( '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 drprft
418 
419 
420 
421  subroutine drprfb(side, trans, direct, storev, m, n, k, v1, v2, ldv2, t, ldt, c, ldc, work, ldwork)
423  implicit none
424 
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, *)
428 
429  integer :: i, j
430  real(kind(1.d0)), parameter :: one=1.d0
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 dcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
459 #if defined(cprec) || defined(zprec)
460  call dlacgv(n, work(1,j), 1)
461 #endif
462  end do
463  ! w := w * v1
464  call drpmm( 'right', 'lower', notran, 'unit', n, k, one, v1, work, ldwork )
465  if( m.gt.k ) then
466  ! w := w + c2'*v2
467  call dgemm( 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 dtrmm( '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 dgemm( 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 drpmm( '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 ) - (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 drprfb
502 
503 
504 
505 
506 
507  subroutine drpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
509  implicit none
510  character :: uplo, trans, diag, side
511  integer :: m, n, ldx
512  real(kind(1.d0)) :: 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 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)
531 
532  else
533 
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)
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 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)
553 
554  else
555 
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)
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 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)
579 
580  else
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)
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 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)
600 
601  else
602 
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)
606 
607  end if
608 
609  end if
610 
611  end if
612 
613 
614 
615  return
616 
617  end subroutine drpmm
618 
619 
620 
621  subroutine drpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
623  implicit none
624  character :: uplo, trans, diag, side
625  integer :: m, n, ldx
626  real(kind(1.d0)) :: 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 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)
645 
646  else
647 
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)
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 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)
667 
668  else
669 
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)
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 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)
693 
694  else
695 
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)
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 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)
715 
716  else
717 
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)
721 
722  end if
723 
724  end if
725 
726  end if
727 
728  return
729 
730  end subroutine drpsm
731 
732 
733  subroutine dqrm_print_rfpf(uplo, n, b)
735  character :: uplo
736  integer :: n
737  real(kind(1.d0)) :: 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 dqrm_print_rfpf
758 
759 
760 
761 end module dqrm_rfpf_mod
subroutine drpsv(uplo, trans, diag, n, a, x)
character, parameter tran
subroutine drpsm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)
This module contains all the error management routines and data.
subroutine drpmv(uplo, trans, diag, n, a, x)
This module contains an implementation of some operations on triangular/trapezoidal matrices stored i...
subroutine dqrm_print_rfpf(uplo, n, b)
subroutine dqrm_to_rfpf(uplo, unit, n, a, lda, b)
character, parameter notran
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)
This module contains various string handling routines.
subroutine drpmm(side, uplo, trans, diag, m, n, alpha, a, x, ldx)