MAGMA  magma-1.4.0 Matrix Algebra on GPU and Multicore Architectures
clarfb_gpu.cpp
Go to the documentation of this file.
1 /*
2  -- MAGMA (version 1.4.0) --
3  Univ. of Tennessee, Knoxville
4  Univ. of California, Berkeley
6  August 2013
7
8  @author Stan Tomov
9  @author Mark Gates
10  @generated c Wed Aug 14 12:16:09 2013
11 */
12 #include "common_magma.h"
13
14 extern "C" magma_int_t
15 magma_clarfb_gpu( char side, char trans, char direct, char storev,
17  const magmaFloatComplex *dV, magma_int_t ldv,
18  const magmaFloatComplex *dT, magma_int_t ldt,
19  magmaFloatComplex *dC, magma_int_t ldc,
20  magmaFloatComplex *dwork, magma_int_t ldwork )
21 {
22 /* -- MAGMA (version 1.4.0) --
23  Univ. of Tennessee, Univ. of California Berkeley
24  August 2013
25
26  Purpose
27  =======
28  CLARFB applies a complex block reflector H or its transpose H^H to a
29  COMPLEX m by n matrix C, from the left.
30
31  Arguments
32  =========
33  SIDE (input) CHARACTER
34  = 'L': apply H or H^H from the Left
35  = 'R': apply H or H^H from the Right
36
37  TRANS (input) CHARACTER
38  = 'N': apply H (No transpose)
39  = 'C': apply H^H (Conjugate transpose)
40
41  DIRECT (input) CHARACTER
42  Indicates how H is formed from a product of elementary
43  reflectors
44  = 'F': H = H(1) H(2) . . . H(k) (Forward)
45  = 'B': H = H(k) . . . H(2) H(1) (Backward)
46
47  STOREV (input) CHARACTER
48  Indicates how the vectors which define the elementary
49  reflectors are stored:
50  = 'C': Columnwise
51  = 'R': Rowwise
52
53  M (input) INTEGER
54  The number of rows of the matrix C.
55
56  N (input) INTEGER
57  The number of columns of the matrix C.
58
59  K (input) INTEGER
60  The order of the matrix T (= the number of elementary
61  reflectors whose product defines the block reflector).
62
63  DV (input) COMPLEX array on the GPU, dimension
64  (LDV,K) if STOREV = 'C'
65  (LDV,M) if STOREV = 'R' and SIDE = 'L'
66  (LDV,N) if STOREV = 'R' and SIDE = 'R'
67  The matrix V. See further details.
68
69  LDV (input) INTEGER
70  The leading dimension of the array V.
71  If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
72  if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
73  if STOREV = 'R', LDV >= K.
74
75  DT (input) COMPLEX array on the GPU, dimension (LDT,K)
76  The triangular k by k matrix T in the representation of the
77  block reflector.
78
79  LDT (input) INTEGER
80  The leading dimension of the array T. LDT >= K.
81
82  DC (input/output) COMPLEX array on the GPU, dimension (LDC,N)
83  On entry, the m by n matrix C.
84  On exit, C is overwritten by H*C, or H^H*C, or C*H, or C*H^H.
85
86  LDC (input) INTEGER
87  The leading dimension of the array C. LDA >= max(1,M).
88
89  WORK (workspace) COMPLEX array, dimension (LDWORK,K)
90
91  LDWORK (input) INTEGER
92  The leading dimension of the array WORK.
93  If SIDE == 'L', LDWORK >= max(1,N);
94  if SIDE == 'R', LDWORK >= max(1,M);
95
96  Further Details
97  ===============
98  The shape of the matrix V and the storage of the vectors which define
99  the H(i) is best illustrated by the following example with n = 5 and
100  k = 3.
101  All elements including 0's and 1's are stored, unlike LAPACK.
102
103  DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
104
105  V = ( 1 0 0 ) V = ( 1 v1 v1 v1 v1 )
106  ( v1 1 0 ) ( 0 1 v2 v2 v2 )
107  ( v1 v2 1 ) ( 0 0 1 v3 v3 )
108  ( v1 v2 v3 )
109  ( v1 v2 v3 )
110
111  DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
112
113  V = ( v1 v2 v3 ) V = ( v1 v1 1 0 0 )
114  ( v1 v2 v3 ) ( v2 v2 v2 1 0 )
115  ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
116  ( 0 1 v3 )
117  ( 0 0 1 )
118
119  =================================================================== */
120
121  magmaFloatComplex c_zero = MAGMA_C_ZERO;
122  magmaFloatComplex c_one = MAGMA_C_ONE;
123  magmaFloatComplex c_neg_one = MAGMA_C_NEG_ONE;
124
125  /* Check input arguments */
126  magma_int_t info = 0;
127  if (m < 0) {
128  info = -5;
129  } else if (n < 0) {
130  info = -6;
131  } else if (k < 0) {
132  info = -7;
133  } else if ( ((storev == 'C' || storev == 'c') && (side == 'L' || side == 'l') && ldv < max(1,m)) ||
134  ((storev == 'C' || storev == 'c') && (side == 'R' || side == 'r') && ldv < max(1,n)) ||
135  ((storev == 'R' || storev == 'r') && ldv < k) ) {
136  info = -9;
137  } else if (ldt < k) {
138  info = -11;
139  } else if (ldc < max(1,m)) {
140  info = -13;
141  } else if ( ((side == 'L' || side == 'l') && ldwork < max(1,n)) ||
142  ((side == 'R' || side == 'r') && ldwork < max(1,m)) ) {
143  info = -15;
144  }
145  if (info != 0) {
146  magma_xerbla( __func__, -(info) );
147  return info;
148  }
149
150  /* Function Body */
151  if (m <= 0 || n <= 0) {
152  return MAGMA_SUCCESS;
153  }
154
155  // opposite of trans
156  char transt;
157  if (trans == 'N' || trans == 'n')
158  transt = MagmaConjTrans;
159  else
160  transt = MagmaNoTrans;
161
162  // whether T is upper or lower triangular
163  char uplo;
164  if (direct == 'F' || direct == 'f')
165  uplo = MagmaUpper;
166  else
167  uplo = MagmaLower;
168
169  // whether V is stored transposed or not
170  char notransV, transV;
171  if (storev == 'C' || storev == 'c') {
172  notransV = MagmaNoTrans;
173  transV = MagmaConjTrans;
174  }
175  else {
176  notransV = MagmaConjTrans;
177  transV = MagmaNoTrans;
178  }
179
180  if ( side == 'l' || side == 'L' ) {
181  // Form H C or H^H C
182  // Comments assume H C. When forming H^H C, T gets transposed via transt.
183
184  // W = C^H V
185  magma_cgemm( MagmaConjTrans, notransV,
186  n, k, m,
187  c_one, dC, ldc,
188  dV, ldv,
189  c_zero, dwork, ldwork);
190
191  // W = W T^H = C^H V T^H
192  magma_ctrmm( MagmaRight, uplo, transt, MagmaNonUnit,
193  n, k,
194  c_one, dT, ldt,
195  dwork, ldwork);
196
197  // C = C - V W^H = C - V T V^H C = (I - V T V^H) C = H C
198  magma_cgemm( notransV, MagmaConjTrans,
199  m, n, k,
200  c_neg_one, dV, ldv,
201  dwork, ldwork,
202  c_one, dC, ldc);
203  }
204  else {
205  // Form C H or C H^H
206  // Comments assume C H. When forming C H^H, T gets transposed via trans.
207
208  // W = C V
209  magma_cgemm( MagmaNoTrans, notransV,
210  m, k, n,
211  c_one, dC, ldc,
212  dV, ldv,
213  c_zero, dwork, ldwork);
214
215  // W = W T = C V T
216  magma_ctrmm( MagmaRight, uplo, trans, MagmaNonUnit,
217  m, k,
218  c_one, dT, ldt,
219  dwork, ldwork);
220
221  // C = C - W V^H = C - C V T V^H = C (I - V T V^H) = C H
222  magma_cgemm( MagmaNoTrans, transV,
223  m, n, k,
224  c_neg_one, dwork, ldwork,
225  dV, ldv,
226  c_one, dC, ldc);
227  }
228
229  return MAGMA_SUCCESS;
230 } /* magma_clarfb */
magma_int_t magma_clarfb_gpu(char side, char trans, char direct, char storev, magma_int_t m, magma_int_t n, magma_int_t k, const magmaFloatComplex *dv, magma_int_t ldv, const magmaFloatComplex *dt, magma_int_t ldt, magmaFloatComplex *dc, magma_int_t ldc, magmaFloatComplex *dwork, magma_int_t ldwork)
Definition: clarfb_gpu.cpp:15
#define __func__
Definition: common_magma.h:65
#define MagmaUpper
Definition: magma.h:61
#define MAGMA_C_NEG_ONE
Definition: magma.h:156
int magma_int_t
Definition: magmablas.h:12
void magma_cgemm(magma_trans_t transA, magma_trans_t transB, magma_int_t m, magma_int_t n, magma_int_t k, magmaFloatComplex alpha, magmaFloatComplex_const_ptr dA, magma_int_t ldda, magmaFloatComplex_const_ptr dB, magma_int_t lddb, magmaFloatComplex beta, magmaFloatComplex_ptr dC, magma_int_t lddc)
#define dwork(dev, i, j)
#define MagmaLower
Definition: magma.h:62
#define dV(m)
void magma_xerbla(const char *srname, magma_int_t info)
Definition: xerbla.cpp:8
#define MagmaConjTrans
Definition: magma.h:59
#define MAGMA_C_ONE
Definition: magma.h:154
#define MagmaNonUnit
Definition: magma.h:65
#define MAGMA_SUCCESS
Definition: magma.h:106
#define dC(dev, i, j)
#define MagmaRight
Definition: magma.h:69
#define dT(m)
void magma_ctrmm(magma_side_t side, magma_uplo_t uplo, magma_trans_t trans, magma_diag_t diag, magma_int_t m, magma_int_t n, magmaFloatComplex alpha, magmaFloatComplex_const_ptr dA, magma_int_t ldda, magmaFloatComplex_ptr dB, magma_int_t lddb)
#define MagmaNoTrans
Definition: magma.h:57
#define max(a, b)
Definition: common_magma.h:82
#define MAGMA_C_ZERO
Definition: magma.h:153