1 /*
2 Copyright (c) 2009-2014, Jack Poulson
3 All rights reserved.
4
5 This file is part of Elemental and is under the BSD 2-Clause License,
6 which can be found in the LICENSE file in the root directory, or at
7 http://opensource.org/licenses/BSD-2-Clause
8 */
9 #pragma once
10 #ifndef ELEM_MULTISHIFTQUASITRSM_HPP
11 #define ELEM_MULTISHIFTQUASITRSM_HPP
12
13 namespace elem {
14
15 template<typename F,Dist shiftColDist,
16 Dist XColDist,Dist XRowDist>
17 inline void
LocalMultiShiftQuasiTrsm(LeftOrRight side,UpperOrLower uplo,Orientation orientation,F alpha,const DistMatrix<F,STAR,STAR> & A,const DistMatrix<F,shiftColDist,STAR> & shifts,DistMatrix<F,XColDist,XRowDist> & X)18 LocalMultiShiftQuasiTrsm
19 ( LeftOrRight side, UpperOrLower uplo, Orientation orientation,
20 F alpha, const DistMatrix<F,STAR,STAR>& A,
21 const DistMatrix<F,shiftColDist,STAR >& shifts,
22 DistMatrix<F, XColDist,XRowDist>& X )
23 {
24 DEBUG_ONLY(
25 CallStackEntry cse("LocalMultiShiftQuasiTrsm");
26 if( (side == LEFT && ( XColDist != STAR ||
27 shiftColDist != XRowDist) ) ||
28 (side == RIGHT && ( XRowDist != STAR ||
29 shiftColDist != XColDist) ) )
30 LogicError
31 ("Dist of RHS and shifts must conform with that of triangle");
32 )
33 // NOTE: Is this prototype available yet?!?
34 MultiShiftQuasiTrsm
35 ( side, uplo, orientation,
36 alpha, A.LockedMatrix(), shifts.LockedMatrix(), X.Matrix() );
37 }
38
39 template<typename Real,Dist shiftColDist,
40 Dist XColDist,Dist XRowDist>
41 inline void
LocalMultiShiftQuasiTrsm(LeftOrRight side,UpperOrLower uplo,Orientation orientation,Complex<Real> alpha,const DistMatrix<Real,STAR,STAR> & A,const DistMatrix<Complex<Real>,shiftColDist,STAR> & shifts,DistMatrix<Real,XColDist,XRowDist> & XReal,DistMatrix<Real,XColDist,XRowDist> & XImag)42 LocalMultiShiftQuasiTrsm
43 ( LeftOrRight side, UpperOrLower uplo, Orientation orientation,
44 Complex<Real> alpha,
45 const DistMatrix<Real, STAR, STAR >& A,
46 const DistMatrix<Complex<Real>,shiftColDist,STAR >& shifts,
47 DistMatrix<Real, XColDist,XRowDist>& XReal,
48 DistMatrix<Real, XColDist,XRowDist>& XImag )
49 {
50 DEBUG_ONLY(
51 CallStackEntry cse("LocalMultiShiftQuasiTrsm");
52 if( (side == LEFT && ( XColDist != STAR ||
53 shiftColDist != XRowDist) ) ||
54 (side == RIGHT && ( XRowDist != STAR ||
55 shiftColDist != XColDist) ) )
56 LogicError
57 ("Dist of RHS and shifts must conform with that of triangle");
58 )
59 // NOTE: Is this prototype available yet?!?
60 MultiShiftQuasiTrsm
61 ( side, uplo, orientation,
62 alpha, A.LockedMatrix(), shifts.LockedMatrix(),
63 XReal.Matrix(), XImag.Matrix() );
64 }
65
66 } // namespace elem
67
68 #include "./MultiShiftQuasiTrsm/LLN.hpp"
69 #include "./MultiShiftQuasiTrsm/LLT.hpp"
70 #include "./MultiShiftQuasiTrsm/LUN.hpp"
71 #include "./MultiShiftQuasiTrsm/LUT.hpp"
72 //#include "./MultiShiftQuasiTrsm/RLN.hpp"
73 //#include "./MultiShiftQuasiTrsm/RLT.hpp"
74 //#include "./MultiShiftQuasiTrsm/RUN.hpp"
75 //#include "./MultiShiftQuasiTrsm/RUT.hpp"
76
77 namespace elem {
78
79 template<typename F>
80 inline void
MultiShiftQuasiTrsm(LeftOrRight side,UpperOrLower uplo,Orientation orientation,F alpha,const Matrix<F> & A,const Matrix<F> & shifts,Matrix<F> & B)81 MultiShiftQuasiTrsm
82 ( LeftOrRight side, UpperOrLower uplo, Orientation orientation,
83 F alpha, const Matrix<F>& A, const Matrix<F>& shifts, Matrix<F>& B )
84 {
85 DEBUG_ONLY(
86 CallStackEntry cse("MultiShiftQuasiTrsm");
87 if( A.Height() != A.Width() )
88 LogicError("A must be square");
89 if( side == LEFT )
90 {
91 if( A.Height() != B.Height() )
92 LogicError("Nonconformal");
93 }
94 else
95 {
96 if( A.Height() != B.Width() )
97 LogicError("Nonconformal");
98 }
99 )
100 Scale( alpha, B );
101 // TODO: Call the single right-hand side algorithm if appropriate
102
103 if( side == LEFT && uplo == LOWER )
104 {
105 if( orientation == NORMAL )
106 msquasitrsm::LLN( A, shifts, B );
107 else
108 msquasitrsm::LLT( orientation, A, shifts, B );
109 }
110 else if( side == LEFT && uplo == UPPER )
111 {
112 if( orientation == NORMAL )
113 msquasitrsm::LUN( A, shifts, B );
114 else
115 msquasitrsm::LUT( orientation, A, shifts, B );
116 }
117 else if( side == RIGHT && uplo == LOWER )
118 {
119 if( orientation == NORMAL )
120 //msquasitrsm::RLN( A, B );
121 LogicError("This case not yet handled");
122 else
123 //msquasitrsm::RLT( orientation, A, B );
124 LogicError("This case not yet handled");
125 }
126 else if( side == RIGHT && uplo == UPPER )
127 {
128 if( orientation == NORMAL )
129 //msquasitrsm::RUN( A, B );
130 LogicError("This case not yet handled");
131 else
132 //msquasitrsm::RUT( orientation, A, B );
133 LogicError("This case not yet handled");
134 }
135 }
136
137 template<typename Real>
138 inline void
MultiShiftQuasiTrsm(LeftOrRight side,UpperOrLower uplo,Orientation orientation,Complex<Real> alpha,const Matrix<Real> & A,const Matrix<Complex<Real>> & shifts,Matrix<Real> & BReal,Matrix<Real> & BImag)139 MultiShiftQuasiTrsm
140 ( LeftOrRight side, UpperOrLower uplo, Orientation orientation,
141 Complex<Real> alpha,
142 const Matrix<Real>& A,
143 const Matrix<Complex<Real>>& shifts,
144 Matrix<Real>& BReal, Matrix<Real>& BImag )
145 {
146 DEBUG_ONLY(
147 CallStackEntry cse("MultiShiftQuasiTrsm");
148 if( A.Height() != A.Width() )
149 LogicError("A must be square");
150 if( side == LEFT )
151 {
152 if( A.Height() != BReal.Height() )
153 LogicError("Nonconformal");
154 }
155 else
156 {
157 if( A.Height() != BReal.Width() )
158 LogicError("Nonconformal");
159 }
160 )
161 Scale( alpha, BReal, BImag );
162 // TODO: Call the single right-hand side algorithm if appropriate
163
164 if( side == LEFT && uplo == LOWER )
165 {
166 if( orientation == NORMAL )
167 //msquasitrsm::LLN( A, shifts, BReal, BImag );
168 LogicError("This case not yet handled");
169 else
170 //msquasitrsm::LLT( orientation, A, shifts, BReal, BImag );
171 LogicError("This case not yet handled");
172 }
173 else if( side == LEFT && uplo == UPPER )
174 {
175 if( orientation == NORMAL )
176 msquasitrsm::LUN( A, shifts, BReal, BImag );
177 else
178 msquasitrsm::LUT( orientation, A, shifts, BReal, BImag );
179 }
180 else if( side == RIGHT && uplo == LOWER )
181 {
182 if( orientation == NORMAL )
183 //msquasitrsm::RLN( A, BReal, BImag );
184 LogicError("This case not yet handled");
185 else
186 //msquasitrsm::RLT( orientation, A, BReal, BImag );
187 LogicError("This case not yet handled");
188 }
189 else if( side == RIGHT && uplo == UPPER )
190 {
191 if( orientation == NORMAL )
192 //msquasitrsm::RUN( A, BReal, BImag );
193 LogicError("This case not yet handled");
194 else
195 //msquasitrsm::RUT( orientation, A, BReal, BImag );
196 LogicError("This case not yet handled");
197 }
198 }
199
200 template<typename F>
201 inline void
MultiShiftQuasiTrsm(LeftOrRight side,UpperOrLower uplo,Orientation orientation,F alpha,const DistMatrix<F> & A,const DistMatrix<F,VR,STAR> & shifts,DistMatrix<F> & B)202 MultiShiftQuasiTrsm
203 ( LeftOrRight side, UpperOrLower uplo, Orientation orientation,
204 F alpha, const DistMatrix<F>& A, const DistMatrix<F,VR,STAR>& shifts,
205 DistMatrix<F>& B )
206 {
207 DEBUG_ONLY(
208 CallStackEntry cse("MultiShiftQuasiTrsm");
209 if( A.Grid() != B.Grid() )
210 LogicError("A and B must use the same grid");
211 if( A.Height() != A.Width() )
212 LogicError("A must be square");
213 if( side == LEFT )
214 {
215 if( A.Height() != B.Height() )
216 LogicError("Nonconformal");
217 }
218 else
219 {
220 if( A.Height() != B.Width() )
221 LogicError("Nonconformal");
222 }
223 )
224 Scale( alpha, B );
225 // TODO: Call the single right-hand side algorithm if appropriate
226
227 //const Int p = B.Grid().Size();
228 if( side == LEFT && uplo == LOWER )
229 {
230 if( orientation == NORMAL )
231 {
232 msquasitrsm::LLNLarge( A, shifts, B );
233 /*
234 if( B.Width() > 5*p )
235 msquasitrsm::LLNLarge( A, shifts, B );
236 else
237 msquasitrsm::LLNMedium( A, shifts, B );
238 */
239 }
240 else
241 {
242 msquasitrsm::LLTLarge( orientation, A, shifts, B );
243 /*
244 if( B.Width() > 5*p )
245 msquasitrsm::LLTLarge( orientation, A, shifts, B );
246 else
247 msquasitrsm::LLTMedium( orientation, A, shifts, B );
248 */
249 }
250 }
251 else if( side == LEFT && uplo == UPPER )
252 {
253 if( orientation == NORMAL )
254 {
255 msquasitrsm::LUNLarge( A, shifts, B );
256 /*
257 if( B.Width() > 5*p )
258 msquasitrsm::LUNLarge( A, shifts, B );
259 else
260 msquasitrsm::LUNMedium( A, shifts, B );
261 */
262 }
263 else
264 {
265 msquasitrsm::LUTLarge( orientation, A, shifts, B );
266 /*
267 if( B.Width() > 5*p )
268 msquasitrsm::LUTLarge( orientation, A, shifts, B );
269 else
270 msquasitrsm::LUTMedium( orientation, A, shifts, B );
271 */
272 }
273 }
274 else if( side == RIGHT && uplo == LOWER )
275 {
276 if( orientation == NORMAL )
277 //msquasitrsm::RLN( A, shifts, B );
278 LogicError("This case not yet handled");
279 else
280 //msquasitrsm::RLT( orientation, A, shifts, B );
281 LogicError("This case not yet handled");
282 }
283 else if( side == RIGHT && uplo == UPPER )
284 {
285 if( orientation == NORMAL )
286 //msquasitrsm::RUN( A, shifts, B );
287 LogicError("This case not yet handled");
288 else
289 //msquasitrsm::RUT( orientation, A, shifts, B );
290 LogicError("This case not yet handled");
291 }
292 }
293
294 template<typename Real>
295 inline void
MultiShiftQuasiTrsm(LeftOrRight side,UpperOrLower uplo,Orientation orientation,Complex<Real> alpha,const DistMatrix<Real> & A,const DistMatrix<Complex<Real>,VR,STAR> & shifts,DistMatrix<Real> & BReal,DistMatrix<Real> & BImag)296 MultiShiftQuasiTrsm
297 ( LeftOrRight side, UpperOrLower uplo, Orientation orientation,
298 Complex<Real> alpha,
299 const DistMatrix<Real>& A,
300 const DistMatrix<Complex<Real>,VR,STAR>& shifts,
301 DistMatrix<Real>& BReal, DistMatrix<Real>& BImag )
302 {
303 DEBUG_ONLY(
304 CallStackEntry cse("MultiShiftQuasiTrsm");
305 if( A.Grid() != BReal.Grid() || BReal.Grid() != BImag.Grid() )
306 LogicError("A and B must use the same grid");
307 if( A.Height() != A.Width() )
308 LogicError("A must be square");
309 if( BReal.Height() != BImag.Height() ||
310 BReal.Width() != BImag.Width() )
311 LogicError("BReal and BImag must be the same size");
312 if( side == LEFT )
313 {
314 if( A.Height() != BReal.Height() )
315 LogicError("Nonconformal");
316 }
317 else
318 {
319 if( A.Height() != BReal.Width() )
320 LogicError("Nonconformal");
321 }
322 )
323 Scale( alpha, BReal, BImag );
324 // TODO: Call the single right-hand side algorithm if appropriate
325
326 const Int p = BReal.Grid().Size();
327 if( side == LEFT && uplo == LOWER )
328 {
329 if( orientation == NORMAL )
330 {
331 if( BReal.Width() > 5*p )
332 //msquasitrsm::LLNLarge( A, shifts, BReal, BImag );
333 LogicError("This case not yet handled");
334 else
335 //msquasitrsm::LLNMedium( A, shifts, BReal, BImag );
336 LogicError("This case not yet handled");
337 }
338 else
339 {
340 if( BReal.Width() > 5*p )
341 //msquasitrsm::LLTLarge( orientation, A, shifts, BReal, BImag );
342 LogicError("This case not yet handled");
343 else
344 //msquasitrsm::LLTMedium( orientation, A, shifts, BReal, BImag );
345 LogicError("This case not yet handled");
346 }
347 }
348 else if( side == LEFT && uplo == UPPER )
349 {
350 if( orientation == NORMAL )
351 {
352 msquasitrsm::LUNLarge( A, shifts, BReal, BImag );
353 /*
354 if( BReal.Width() > 5*p )
355 msquasitrsm::LUNLarge( A, shifts, BReal, BImag );
356 else
357 msquasitrsm::LUNMedium( A, shifts, BReal, BImag );
358 */
359 }
360 else
361 {
362 msquasitrsm::LUTLarge( orientation, A, shifts, BReal, BImag );
363 /*
364 if( BReal.Width() > 5*p )
365 msquasitrsm::LUTLarge( orientation, A, shifts, BReal, BImag );
366 else
367 msquasitrsm::LUTMedium( orientation, A, shifts, BReal, BImag );
368 */
369 }
370 }
371 else if( side == RIGHT && uplo == LOWER )
372 {
373 if( orientation == NORMAL )
374 //msquasitrsm::RLN( A, shifts, BReal, BImag );
375 LogicError("This case not yet handled");
376 else
377 //msquasitrsm::RLT( orientation, A, shifts, BReal, BImag );
378 LogicError("This case not yet handled");
379 }
380 else if( side == RIGHT && uplo == UPPER )
381 {
382 if( orientation == NORMAL )
383 //msquasitrsm::RUN( A, shifts, BReal, BImag );
384 LogicError("This case not yet handled");
385 else
386 //msquasitrsm::RUT( orientation, A, shifts, BReal, BImag );
387 LogicError("This case not yet handled");
388 }
389 }
390
391 } // namespace elem
392
393 #endif // ifndef ELEM_MULTISHIFTQUASITRSM_HPP
394