QR_MUMPS
 All Classes Files Functions Variables Enumerations Enumerator Pages
dqrm_matmul.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 
51 subroutine dqrm_matmul2d(qrm_mat, transp, alpha, x, beta, y)
52 
53  use dqrm_spmat_mod
54  use qrm_string_mod
55  implicit none
56 
57  type(dqrm_spmat_type) :: qrm_mat
58  real(kind(1.d0)), intent(out) :: y(:,:)
59  real(kind(1.d0)), intent(in) :: x(:,:)
60  real(kind(1.d0)), intent(in) :: alpha, beta
61  character(len=*) :: transp
62 
63  integer :: nb, nrhs, j, c, r, k, i, rhs_nthreads
64 
65  call qrm_get(qrm_mat, 'qrm_rhsnb', nb)
66  call qrm_get(qrm_mat, 'qrm_rhsnthreads', rhs_nthreads)
67  nrhs = size(x,2)
68  if(nb.le.0) nb = nrhs
69 
70  ! TODO: add checks for sizes etc.
71 
72  if(beta .eq. 0.d0) then
73  y = 0.d0
74  else
75  y = beta*y
76  end if
77 
78  ! shortcut
79  if(alpha .eq. 0.d0) then
80  return
81  end if
82 
83  !$omp parallel do num_threads(rhs_nthreads) private(i, k, r, c)
84  do k=1, nrhs, nb
85  do i=1, qrm_mat%nz
86  if((qrm_str_tolower(transp(1:1)) .eq. 't') .or. (qrm_str_tolower(transp(1:1)) .eq. 'c')) then
87  c = qrm_mat%irn(i)
88  r = qrm_mat%jcn(i)
89  y(r,k:min(nrhs,k+nb-1)) = y(r,k:min(nrhs,k+nb-1))+alpha*(qrm_mat%val(i))*x(c,k:min(nrhs,k+nb-1))
90  else
91  r = qrm_mat%irn(i)
92  c = qrm_mat%jcn(i)
93  y(r,k:min(nrhs,k+nb-1)) = y(r,k:min(nrhs,k+nb-1))+alpha*qrm_mat%val(i)*x(c,k:min(nrhs,k+nb-1))
94  end if
95  end do
96  end do
97  !$omp end parallel do
98 
99  return
100 
101 end subroutine dqrm_matmul2d
102 
103 
104 
122 subroutine dqrm_matmul1d(qrm_mat, transp, alpha, x, beta, y)
123 
124  use dqrm_spmat_mod
125  use qrm_string_mod
126  implicit none
127 
128  type(dqrm_spmat_type) :: qrm_mat
129  real(kind(1.d0)), intent(out) :: y(:)
130  real(kind(1.d0)), intent(in) :: x(:)
131  real(kind(1.d0)), intent(in) :: alpha, beta
132  character(len=*) :: transp
133 
134  integer :: c, r, i, n
135 
136  n = size(x,1)
137 
138  ! TODO: add checks for sizes etc.
139 
140  if(beta .eq. 0.d0) then
141  y = 0.d0
142  else
143  y = beta*y
144  end if
145 
146  ! shortcut
147  if(alpha .eq. 0.d0) then
148  return
149  end if
150 
151  do i=1, qrm_mat%nz
152  if((qrm_str_tolower(transp(1:1)) .eq. 't') .or. (qrm_str_tolower(transp(1:1)) .eq. 'c')) then
153  c = qrm_mat%irn(i)
154  r = qrm_mat%jcn(i)
155  y(r) = y(r)+alpha*(qrm_mat%val(i))*x(c)
156  else
157  r = qrm_mat%irn(i)
158  c = qrm_mat%jcn(i)
159  y(r) = y(r)+alpha*qrm_mat%val(i)*x(c)
160  end if
161  end do
162 
163  return
164 
165 end subroutine dqrm_matmul1d
Generif interface for the ::dqrm_pgeti, ::dqrm_pgetr and.
This module contains the definition of the basic sparse matrix type and of the associated methods...
subroutine dqrm_matmul2d(qrm_mat, transp, alpha, x, beta, y)
This subroutine computes the product y=beta*y+alpha*op(A)*x where op(A) is either A or A' depending o...
Definition: dqrm_matmul.F90:51
This type defines the data structure used to store a matrix.
subroutine dqrm_matmul1d(qrm_mat, transp, alpha, x, beta, y)
This subroutine computes the product y=beta*y+alpha*op(A)*x where op(A) is either A or A' depending o...
This module contains various string handling routines.
int i
Definition: secs.c:40