QR_MUMPS
 All Classes Files Functions Variables Enumerations Enumerator Pages
qrm_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 _qrm_solve2d(qrm_mat, transp, b, x)
51 
52  use _qrm_spmat_mod
53  use qrm_error_mod
54  use _qrm_fdata_mod
55  use qrm_string_mod
56  use _qrm_utils_mod
57  use _qrm_solve_mod, protect => _qrm_solve2d
58  implicit none
59 
60  type(_qrm_spmat_type), target :: qrm_mat
61  _qrm_data, intent(inout) :: b(:,:)
62  _qrm_data, 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 = _qrm_zero
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 _qrm_solve2d
114 
115 
116 
128 subroutine _qrm_solve1d(qrm_mat, transp, b, x)
129 
130  use _qrm_spmat_mod
131  use qrm_error_mod
132  use _qrm_fdata_mod
133  use qrm_string_mod
134  use _qrm_utils_mod
135  use _qrm_solve_mod, protect => _qrm_solve1d
136  use qrm_error_mod
137  implicit none
138 
139  type(_qrm_spmat_type), target :: qrm_mat
140  _qrm_data, intent(in) :: b(:)
141  _qrm_data, intent(out) :: x(:)
142  character(len=*) :: transp
143 
144  _qrm_data, 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 = _qrm_zero
159 
160  n = size(b,1)
161  call _qrm_remap_pnt(b, pntb, n)
162  n = size(x,1)
163  call _qrm_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 _qrm_solve1d
184 
This module contains all the error management routines and data.
This module contains the definition of all the data related to the factorization phase.
subroutine _qrm_remap_pnt(arr1d, pnt2d, n)
This function makes a 2D pointer point to a 1D array.
subroutine qrm_err_act_restore(err_act)
Restores the value of the qrm_err_act variable.
This module contains generic interfaces for a number of auxiliary tools.
subroutine qrm_err_check()
This subroutine checks the errors stack. If something is found all the entries in the stack are poppe...
subroutine _qrm_solve1d(qrm_mat, transp, b, x)
This function solves for R or R' against a single vector.
Definition: qrm_solve.F90:128
Generic interface for the ::_qrm_solve_r routine.
Generic interface for the ::_qrm_solve_rt routine.
This module contains all the interfaces for the typed routines in the solve phase.
This type defines the data structure used to store a matrix.
This module contains the definition of the basic sparse matrix type and of the associated methods...
Generif interface for the ::_qrm_pgeti, ::_qrm_pgetr and.
This module contains various string handling routines.
int i
Definition: secs.c:40
subroutine _qrm_solve2d(qrm_mat, transp, b, x)
This function solves for R or R' against multiple vectors.
Definition: qrm_solve.F90:50
subroutine qrm_err_push(code, sub, ied, aed)
This subroutine pushes an error on top of the stack.
subroutine qrm_err_act_save(err_act)
Saves a copy of the qrm_err_act variable.