QR_MUMPS
dqrm_solve.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 "prec.h"
36 #include "qrm_common.h"
37 
38 
50 subroutine dqrm_solve2d(qrm_mat, transp, b, x)
51 
52  use dqrm_spmat_mod
53  use qrm_error_mod
54  use dqrm_fdata_mod
55  use qrm_string_mod
56  use dqrm_utils_mod
57  use dqrm_solve_mod, protect => dqrm_solve2d
58  implicit none
59 
60  type(dqrm_spmat_type), target :: qrm_mat
61  real(kind(1.d0)), intent(inout) :: b(:,:)
62  real(kind(1.d0)), intent(out) :: x(:,:)
63  character(len=*) :: transp
64 
65  integer :: i, nb, nrhs, rhs_nthreads
66  ! error management
67  integer :: err_act
68  character(len=*), parameter :: name='qrm_solve'
69 
70  call qrm_err_act_save(err_act)
71 
72  if(.not. qrm_mat%fdata%ok) then
73  call qrm_err_push(14, 'qrm_solve')
74  return
75  end if
76 
77  ! FIXME x should not be initialized to 0
78  x = 0.d0
79 
80  ! blocking to deal with multiple rhs
81  call qrm_get(qrm_mat, 'qrm_rhsnb', nb)
82  call qrm_get(qrm_mat, 'qrm_rhsnthreads', rhs_nthreads)
83  nrhs = size(b,2)
84  if(nb.le.0) nb = nrhs
85 
86  !$ call omp_set_nested(.true.)
87  if(qrm_str_tolower(transp(1:1)) .eq. 't') then
88  !$omp parallel do num_threads(rhs_nthreads) private(i)
89  do i=1, nrhs, nb
90  call qrm_solve_rt(qrm_mat, b(:,i:min(nrhs,i+nb-1)), x(:,i:min(nrhs,i+nb-1)))
91  end do
92  !$omp end parallel do
93  __qrm_check_ret(name,'qrm_solve_rt',9999)
94  else
95  !$omp parallel do num_threads(rhs_nthreads) private(i)
96  do i=1, nrhs, nb
97  call qrm_solve_r(qrm_mat, b(:,i:min(nrhs,i+nb-1)), x(:,i:min(nrhs,i+nb-1)))
98  end do
99  !$omp end parallel do
100  __qrm_check_ret(name,'qrm_solve_r',9999)
101  end if
102 
103  call qrm_err_act_restore(err_act)
104  return
105 
106 9999 continue ! error management
107  call qrm_err_act_restore(err_act)
108  if(err_act .eq. qrm_abort_) then
109  call qrm_err_check()
110  end if
111  return
112 
113 end subroutine dqrm_solve2d
114 
115 
116 
128 subroutine dqrm_solve1d(qrm_mat, transp, b, x)
130  use dqrm_spmat_mod
131  use qrm_error_mod
132  use dqrm_fdata_mod
133  use qrm_string_mod
134  use dqrm_utils_mod
135  use dqrm_solve_mod, protect => dqrm_solve1d
136  use qrm_error_mod
137  implicit none
138 
139  type(dqrm_spmat_type), target :: qrm_mat
140  real(kind(1.d0)), intent(in) :: b(:)
141  real(kind(1.d0)), intent(out) :: x(:)
142  character(len=*) :: transp
143 
144  real(kind(1.d0)), pointer :: pntb(:,:), pntx(:,:)
145  integer :: n
146  ! error management
147  integer :: err_act
148  character(len=*), parameter :: name='qrm_solve1d'
149 
150  call qrm_err_act_save(err_act)
151 
152  if( qrm_mat%fdata%done .eq. 0) then
153  call qrm_err_push(14, 'qrm_solve')
154  goto 9999
155  end if
156 
157  ! FIXME x should not be initialized to 0
158  x = 0.d0
159 
160  n = size(b,1)
161  call dqrm_remap_pnt(b, pntb, n)
162  n = size(x,1)
163  call dqrm_remap_pnt(x, pntx, n)
164 
165  if(qrm_str_tolower(transp(1:1)) .eq. 't') then
166  call qrm_solve_rt(qrm_mat, pntb, pntx)
167  __qrm_check_ret(name,'qrm_solve_rt',9999)
168  else
169  call qrm_solve_r(qrm_mat, pntb, pntx)
170  __qrm_check_ret(name,'qrm_solve_r',9999)
171  end if
172 
173  call qrm_err_act_restore(err_act)
174  return
175 
176 9999 continue ! error management
177  call qrm_err_act_restore(err_act)
178  if(err_act .eq. qrm_abort_) then
179  call qrm_err_check()
180  end if
181  return
182 
183 end subroutine dqrm_solve1d
184 
This module contains generic interfaces for a number of auxiliary tools.
subroutine dqrm_solve2d(qrm_mat, transp, b, x)
This function solves for R or R' against multiple vectors.
Definition: dqrm_solve.F90:51
subroutine qrm_err_push(code, sub, ied, aed)
This subroutine pushes an error on top of the stack.
Generif interface for the ::dqrm_pgeti, ::dqrm_pgetr and.
Generic interface for the ::dqrm_solve_r routine.
subroutine qrm_err_act_save(err_act)
Saves a copy of the qrm_err_act variable.
This module contains all the error management routines and data.
This module contains the definition of the basic sparse matrix type and of the associated methods...
integer, parameter qrm_abort_
Possible actions to be performed upon detection of an error.
subroutine qrm_err_check()
This subroutine checks the errors stack. If something is found all the entries in the stack are poppe...
This type defines the data structure used to store a matrix.
This module contains the definition of all the data related to the factorization phase.
subroutine dqrm_remap_pnt(arr1d, pnt2d, n)
This function makes a 2D pointer point to a 1D array.
Generic interface for the ::dqrm_solve_rt routine.
subroutine dqrm_solve1d(qrm_mat, transp, b, x)
This function solves for R or R' against a single vector.
Definition: dqrm_solve.F90:129
This module contains all the interfaces for the typed routines in the solve phase.
This module contains various string handling routines.
subroutine qrm_err_act_restore(err_act)
Restores the value of the qrm_err_act variable.