Lapack_LeastSquares.cxx
1 // Copyright (C) 2003-2009 Marc Duruflé
2 // Copyright (C) 2001-2009 Vivien Mallet
3 //
4 // This file is part of the linear-algebra library Seldon,
5 // http://seldon.sourceforge.net/.
6 //
7 // Seldon is free software; you can redistribute it and/or modify it under the
8 // terms of the GNU Lesser General Public License as published by the Free
9 // Software Foundation; either version 2.1 of the License, or (at your option)
10 // any later version.
11 //
12 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
13 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 // more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public License
18 // along with Seldon. If not, see http://www.gnu.org/licenses/.
19 
20 
21 #ifndef SELDON_FILE_LAPACK_LEAST_SQUARES_CXX
22 
23 
24 #include "Lapack_LeastSquares.hxx"
25 
26 
27 namespace Seldon
28 {
29 
30 
32  // GETQR //
33 
34 
35  /*** ColMajor ***/
36 
37 
38  template<class Prop0, class Allocator0,
39  class Allocator1>
40  void GetQR(Matrix<float, Prop0, ColMajor, Allocator0>& A,
41  Vector<float, VectFull, Allocator1>& tau,
42  LapackInfo& info)
43  {
44  int m = A.GetM();
45  int n = A.GetN();
46  int lwork = max(m,n);
47  Vector<float, VectFull, Allocator1> work(lwork);
48  tau.Reallocate(min(m, n));
49  sgeqrf_(&m, &n, A.GetData(), &m, tau.GetData(),
50  work.GetData(), &lwork, &info.GetInfoRef());
51  }
52 
53 
54  template<class Prop0, class Allocator0,
55  class Allocator1>
56  void GetQR(Matrix<double, Prop0, ColMajor, Allocator0>& A,
57  Vector<double, VectFull, Allocator1>& tau,
58  LapackInfo& info)
59  {
60  int m = A.GetM();
61  int n = A.GetN();
62  int lwork = max(m,n);
63  Vector<double, VectFull, Allocator1> work(lwork);
64  tau.Reallocate(min(m, n));
65  dgeqrf_(&m, &n, A.GetData(), &m, tau.GetData(),
66  work.GetData(), &lwork, &info.GetInfoRef());
67  }
68 
69 
70  template<class Prop0, class Allocator0,
71  class Allocator1>
72  void GetQR(Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
73  Vector<complex<float>, VectFull, Allocator1>& tau,
74  LapackInfo& info)
75  {
76  int m = A.GetM();
77  int n = A.GetN();
78  int lwork = max(m,n);
79  Vector<complex<float>, VectFull, Allocator1> work(lwork);
80  tau.Reallocate(min(m, n));
81  cgeqrf_(&m, &n, A.GetData(), &m, tau.GetData(),
82  work.GetData(), &lwork, &info.GetInfoRef());
83  }
84 
85 
86  template<class Prop0, class Allocator0,
87  class Allocator1>
88  void GetQR(Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
89  Vector<complex<double>, VectFull, Allocator1>& tau,
90  LapackInfo& info)
91  {
92  int m = A.GetM();
93  int n = A.GetN();
94  int lwork = max(m,n);
95  Vector<complex<double>, VectFull, Allocator1> work(lwork);
96  tau.Reallocate(min(m, n));
97  zgeqrf_(&m, &n, A.GetData(), &m, tau.GetData(),
98  work.GetData(), &lwork, &info.GetInfoRef());
99  }
100 
101 
102  /*** RowMajor ***/
103 
104 
105  template<class Prop0, class Allocator0,
106  class Allocator1>
107  void GetQR(Matrix<float, Prop0, RowMajor, Allocator0>& A,
108  Vector<float, VectFull, Allocator1>& tau,
109  LapackInfo& info)
110  {
111  int m = A.GetM();
112  int n = A.GetN();
113  int lwork = max(m,n);
114  Vector<float, VectFull, Allocator1> work(lwork);
115  tau.Reallocate(min(m, n));
116  // Factorization LQ of A^t.
117  sgelqf_(&n, &m, A.GetData(), &n, tau.GetData(),
118  work.GetData(), &lwork, &info.GetInfoRef());
119  }
120 
121 
122  template<class Prop0, class Allocator0,
123  class Allocator1>
124  void GetQR(Matrix<double, Prop0, RowMajor, Allocator0>& A,
125  Vector<double, VectFull, Allocator1>& tau,
126  LapackInfo& info)
127  {
128  int m = A.GetM();
129  int n = A.GetN();
130  int lwork = max(m,n);
131  Vector<double, VectFull, Allocator1> work(lwork);
132  tau.Reallocate(min(m, n));
133  // Factorization LQ of A^t.
134  dgelqf_(&n, &m, A.GetData(), &n, tau.GetData(),
135  work.GetData(), &lwork, &info.GetInfoRef());
136  }
137 
138 
139  template<class Prop0, class Allocator0,
140  class Allocator1>
141  void GetQR(Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
142  Vector<complex<float>, VectFull, Allocator1>& tau,
143  LapackInfo& info)
144  {
145  int m = A.GetM();
146  int n = A.GetN();
147  int lwork = max(m,n);
148  Vector<complex<float>, VectFull, Allocator1> work(lwork);
149  tau.Reallocate(min(m, n));
150  // Factorization LQ of A^t.
151  cgelqf_(&n, &m, A.GetData(), &n, tau.GetData(),
152  work.GetData(), &lwork, &info.GetInfoRef());
153  }
154 
155 
156  template<class Prop0, class Allocator0,
157  class Allocator1>
158  void GetQR(Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
159  Vector<complex<double>, VectFull, Allocator1>& tau,
160  LapackInfo& info)
161  {
162  int m = A.GetM();
163  int n = A.GetN();
164  int lwork = max(m,n);
165  Vector<complex<double>, VectFull, Allocator1> work(lwork);
166  tau.Reallocate(min(m, n));
167  // Factorization LQ of A^t.
168  zgelqf_(&n, &m, A.GetData(), &n, tau.GetData(),
169  work.GetData(), &lwork, &info.GetInfoRef());
170  }
171 
172 
173  // GETQR //
175 
176 
178  // GETQR_PIVOT //
179 
180 
181  /*** ColMajor ***/
182 
183 
184  template<class Prop0, class Allocator0,
185  class Allocator1>
186  void GetQR_Pivot(Matrix<double, Prop0, ColMajor, Allocator0>& A,
187  Vector<double, VectFull, Allocator1>& tau,
188  Vector<int>& ipivot, LapackInfo& info)
189  {
190  int m = A.GetM();
191  int n = A.GetN();
192  int lwork = 4 * max(m, n);
193  ipivot.Fill(0);
194  Vector<double, VectFull, Allocator1> work(lwork);
195  tau.Reallocate(min(m, n));
196  dgeqp3_(&m, &n, A.GetData(), &m, ipivot.GetData(), tau.GetData(),
197  work.GetData(), &lwork, &info.GetInfoRef());
198  }
199 
200 
201  // GETQR_PIVOT //
203 
204 
206  // GETQ_FROMQR //
207 
208 
209  /*** ColMajor ***/
210 
211  template<class Prop0, class Allocator0,
212  class Allocator1>
213  void GetQ_FromQR(Matrix<float, Prop0, ColMajor, Allocator0>& A,
214  Vector<float, VectFull, Allocator1>& tau,
215  LapackInfo& info)
216  {
217  int m = A.GetM();
218  int n = A.GetN();
219  int lwork = 2 * max(m, n);
220  int k = min(m, n);
221 
222 #ifdef SELDON_CHECK_DIMENSIONS
223  if (tau.GetM() < k)
224  throw WrongDim("GetQ_FromQR", "tau not large enough");
225 #endif
226 
227  Vector<float, VectFull, Allocator1> work(lwork);
228  sorgqr_(&m, &k, &k, A.GetData(), &m, tau.GetData(),
229  work.GetData(), &lwork, &info.GetInfoRef());
230 
231  if (m < n)
232  A.Resize(m, m);
233  }
234 
235 
236  template<class Prop0, class Allocator0,
237  class Allocator1>
238  void GetQ_FromQR(Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
239  Vector<complex<float>, VectFull, Allocator1>& tau,
240  LapackInfo& info)
241  {
242  int m = A.GetM();
243  int n = A.GetN();
244  int lwork = 2 * max(m, n);
245  int k = min(m, n);
246 
247 #ifdef SELDON_CHECK_DIMENSIONS
248  if (tau.GetM() < k)
249  throw WrongDim("GetQ_FromQR", "tau not large enough");
250 #endif
251 
252  Vector<complex<float>, VectFull, Allocator1> work(lwork);
253  cungqr_(&m, &k, &k, A.GetDataVoid(), &m, tau.GetData(),
254  work.GetDataVoid(), &lwork, &info.GetInfoRef());
255 
256  if (m < n)
257  A.Resize(m, m);
258  }
259 
260 
261  template<class Prop0, class Allocator0,
262  class Allocator1>
263  void GetQ_FromQR(Matrix<double, Prop0, ColMajor, Allocator0>& A,
264  Vector<double, VectFull, Allocator1>& tau,
265  LapackInfo& info)
266  {
267  int m = A.GetM();
268  int n = A.GetN();
269  int lwork = 2 * max(m, n);
270  int k = min(m, n);
271 
272 #ifdef SELDON_CHECK_DIMENSIONS
273  if (tau.GetM() < k)
274  throw WrongDim("GetQ_FromQR", "tau not large enough");
275 #endif
276 
277  Vector<double, VectFull, Allocator1> work(lwork);
278  dorgqr_(&m, &k, &k, A.GetData(), &m, tau.GetData(),
279  work.GetData(), &lwork, &info.GetInfoRef());
280 
281  if (m < n)
282  A.Resize(m, m);
283  }
284 
285 
286  template<class Prop0, class Allocator0,
287  class Allocator1>
288  void GetQ_FromQR(Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
289  Vector<complex<double>, VectFull, Allocator1>& tau,
290  LapackInfo& info)
291  {
292  int m = A.GetM();
293  int n = A.GetN();
294  int lwork = 2 * max(m, n);
295  int k = min(m, n);
296 
297 #ifdef SELDON_CHECK_DIMENSIONS
298  if (tau.GetM() < k)
299  throw WrongDim("GetQ_FromQR", "tau not large enough");
300 #endif
301 
302  Vector<complex<double>, VectFull, Allocator1> work(lwork);
303  zungqr_(&m, &k, &k, A.GetDataVoid(), &m, tau.GetData(),
304  work.GetDataVoid(), &lwork, &info.GetInfoRef());
305 
306  if (m < n)
307  A.Resize(m, m);
308  }
309 
310 
311  template<class Prop0, class Allocator0,
312  class Allocator1>
313  void GetQ_FromLQ(Matrix<float, Prop0, ColMajor, Allocator0>& A,
314  Vector<float, VectFull, Allocator1>& tau,
315  LapackInfo& info)
316  {
317  int m = A.GetM();
318  int n = A.GetN();
319  int lwork = 2 * max(m, n);
320 
321 #ifdef SELDON_CHECK_DIMENSIONS
322  if (tau.GetM() < min(m, n))
323  throw WrongDim("GetQ_FromLQ", "tau not large enough");
324 #endif
325 
326  int k = min(m, n);
327  Vector<float, VectFull, Allocator1> work(lwork);
328  sorglq_(&k, &n, &k, A.GetData(), &m, tau.GetData(),
329  work.GetData(), &lwork, &info.GetInfoRef());
330 
331  if (m > n)
332  A.Resize(n, n);
333  }
334 
335 
336  template<class Prop0, class Allocator0,
337  class Allocator1>
338  void GetQ_FromLQ(Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
339  Vector<complex<float>, VectFull, Allocator1>& tau,
340  LapackInfo& info)
341  {
342  int m = A.GetM();
343  int n = A.GetN();
344  int lwork = 2 * max(m, n);
345 
346 #ifdef SELDON_CHECK_DIMENSIONS
347  if (tau.GetM() < min(m, n))
348  throw WrongDim("GetQ_FromLQ", "tau not large enough");
349 #endif
350 
351  int k = min(m, n);
352  Vector<complex<float>, VectFull, Allocator1> work(lwork);
353  cunglq_(&k, &n, &k, A.GetDataVoid(), &m, tau.GetDataVoid(),
354  work.GetDataVoid(), &lwork, &info.GetInfoRef());
355 
356  if (m > n)
357  A.Resize(n, n);
358  }
359 
360 
361  template<class Prop0, class Allocator0,
362  class Allocator1>
363  void GetQ_FromLQ(Matrix<double, Prop0, ColMajor, Allocator0>& A,
364  Vector<double, VectFull, Allocator1>& tau,
365  LapackInfo& info)
366  {
367  int m = A.GetM();
368  int n = A.GetN();
369  int lwork = 2 * max(m, n);
370 
371 #ifdef SELDON_CHECK_DIMENSIONS
372  if (tau.GetM() < min(m, n))
373  throw WrongDim("GetQ_FromLQ", "tau not large enough");
374 #endif
375 
376  int k = min(m, n);
377  Vector<double, VectFull, Allocator1> work(lwork);
378  dorglq_(&k, &n, &k, A.GetData(), &m, tau.GetData(),
379  work.GetData(), &lwork, &info.GetInfoRef());
380 
381  if (m > n)
382  A.Resize(n, n);
383  }
384 
385 
386  template<class Prop0, class Allocator0,
387  class Allocator1>
388  void GetQ_FromLQ(Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
389  Vector<complex<double>, VectFull, Allocator1>& tau,
390  LapackInfo& info)
391  {
392  int m = A.GetM();
393  int n = A.GetN();
394  int lwork = 2 * max(m, n);
395 
396 #ifdef SELDON_CHECK_DIMENSIONS
397  if (tau.GetM() < min(m, n))
398  throw WrongDim("GetQ_FromLQ", "tau not large enough");
399 #endif
400 
401  int k = min(m, n);
402  Vector<complex<double>, VectFull, Allocator1> work(lwork);
403  zunglq_(&k, &n, &k, A.GetDataVoid(), &m, tau.GetDataVoid(),
404  work.GetDataVoid(), &lwork, &info.GetInfoRef());
405 
406  if (m > n)
407  A.Resize(n, n);
408  }
409 
410 
411  template<class Prop0, class Allocator0,
412  class Allocator1>
413  void GetQ_FromQR(Matrix<float, Prop0, RowMajor, Allocator0>& A,
414  Vector<float, VectFull, Allocator1>& tau,
415  LapackInfo& info)
416  {
417  int m = A.GetM();
418  int n = A.GetN();
419  int lwork = 2 * max(m, n);
420  int k = min(m, n);
421 
422 #ifdef SELDON_CHECK_DIMENSIONS
423  if (tau.GetM() < k)
424  throw WrongDim("GetQ_FromQR", "tau not large enough");
425 #endif
426 
427  Vector<float, VectFull, Allocator1> work(lwork);
428  sorglq_(&k, &m, &k, A.GetData(), &n, tau.GetData(),
429  work.GetData(), &lwork, &info.GetInfoRef());
430 
431  if (n > m)
432  A.Resize(m, m);
433 
434  }
435 
436 
437  template<class Prop0, class Allocator0,
438  class Allocator1>
439  void GetQ_FromQR(Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
440  Vector<complex<float>, VectFull, Allocator1>& tau,
441  LapackInfo& info)
442  {
443  int m = A.GetM();
444  int n = A.GetN();
445  int lwork = 2 * max(m, n);
446  int k = min(m, n);
447 
448 #ifdef SELDON_CHECK_DIMENSIONS
449  if (tau.GetM() < k)
450  throw WrongDim("GetQ_FromQR", "tau not large enough");
451 #endif
452 
453  Vector<complex<double>, VectFull, Allocator1> work(lwork);
454  cunglq_(&k, &m, &k, A.GetDataVoid(), &n, tau.GetData(),
455  work.GetDataVoid(), &lwork, &info.GetInfoRef());
456 
457  if (n > m)
458  A.Resize(m, m);
459  }
460 
461 
462  template<class Prop0, class Allocator0,
463  class Allocator1>
464  void GetQ_FromQR(Matrix<double, Prop0, RowMajor, Allocator0>& A,
465  Vector<double, VectFull, Allocator1>& tau,
466  LapackInfo& info)
467  {
468  int m = A.GetM();
469  int n = A.GetN();
470  int lwork = 2 * max(m, n);
471  int k = min(m, n);
472 
473 #ifdef SELDON_CHECK_DIMENSIONS
474  if (tau.GetM() < k)
475  throw WrongDim("GetQ_FromQR", "tau not large enough");
476 #endif
477 
478  Vector<double, VectFull, Allocator1> work(lwork);
479  dorglq_(&k, &m, &k, A.GetData(), &n, tau.GetData(),
480  work.GetData(), &lwork, &info.GetInfoRef());
481 
482  if (n > m)
483  A.Resize(m, m);
484 
485  }
486 
487 
488  template<class Prop0, class Allocator0,
489  class Allocator1>
490  void GetQ_FromQR(Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
491  Vector<complex<double>, VectFull, Allocator1>& tau,
492  LapackInfo& info)
493  {
494  int m = A.GetM();
495  int n = A.GetN();
496  int lwork = 2 * max(m, n);
497  int k = min(m, n);
498 
499 #ifdef SELDON_CHECK_DIMENSIONS
500  if (tau.GetM() < k)
501  throw WrongDim("GetQ_FromQR", "tau not large enough");
502 #endif
503 
504  Vector<complex<double>, VectFull, Allocator1> work(lwork);
505  zunglq_(&k, &m, &k, A.GetDataVoid(), &n, tau.GetData(),
506  work.GetDataVoid(), &lwork, &info.GetInfoRef());
507 
508  if (n > m)
509  A.Resize(m, m);
510  }
511 
512 
513  template<class Prop0, class Allocator0,
514  class Allocator1>
515  void GetQ_FromLQ(Matrix<float, Prop0, RowMajor, Allocator0>& A,
516  Vector<float, VectFull, Allocator1>& tau,
517  LapackInfo& info)
518  {
519  int m = A.GetM();
520  int n = A.GetN();
521  int lwork = 2 * max(m, n);
522 
523 #ifdef SELDON_CHECK_DIMENSIONS
524  if (tau.GetM() < min(m, n))
525  throw WrongDim("GetQ_FromLQ", "tau not large enough");
526 #endif
527 
528  int k = min(m, n);
529  Vector<float, VectFull, Allocator1> work(lwork);
530  sorgqr_(&n, &k, &k, A.GetData(), &n, tau.GetData(),
531  work.GetData(), &lwork, &info.GetInfoRef());
532 
533  if (m > n)
534  A.Resize(n, n);
535  }
536 
537 
538  template<class Prop0, class Allocator0,
539  class Allocator1>
540  void GetQ_FromLQ(Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
541  Vector<complex<float>, VectFull, Allocator1>& tau,
542  LapackInfo& info)
543  {
544  int m = A.GetM();
545  int n = A.GetN();
546  int lwork = 2 * max(m, n);
547 
548 #ifdef SELDON_CHECK_DIMENSIONS
549  if (tau.GetM() < min(m, n))
550  throw WrongDim("GetQ_FromLQ", "tau not large enough");
551 #endif
552 
553  int k = min(m, n);
554  Vector<complex<float>, VectFull, Allocator1> work(lwork);
555  cungqr_(&n, &k, &k, A.GetDataVoid(), &n, tau.GetDataVoid(),
556  work.GetDataVoid(), &lwork, &info.GetInfoRef());
557 
558  if (m > n)
559  A.Resize(n, n);
560  }
561 
562 
563  template<class Prop0, class Allocator0,
564  class Allocator1>
565  void GetQ_FromLQ(Matrix<double, Prop0, RowMajor, Allocator0>& A,
566  Vector<double, VectFull, Allocator1>& tau,
567  LapackInfo& info)
568  {
569  int m = A.GetM();
570  int n = A.GetN();
571  int lwork = 2 * max(m, n);
572 
573 #ifdef SELDON_CHECK_DIMENSIONS
574  if (tau.GetM() < min(m, n))
575  throw WrongDim("GetQ_FromLQ", "tau not large enough");
576 #endif
577 
578  int k = min(m, n);
579  Vector<double, VectFull, Allocator1> work(lwork);
580  dorgqr_(&n, &k, &k, A.GetData(), &n, tau.GetData(),
581  work.GetData(), &lwork, &info.GetInfoRef());
582 
583  if (m > n)
584  A.Resize(n, n);
585  }
586 
587 
588  template<class Prop0, class Allocator0,
589  class Allocator1>
590  void GetQ_FromLQ(Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
591  Vector<complex<double>, VectFull, Allocator1>& tau,
592  LapackInfo& info)
593  {
594  int m = A.GetM();
595  int n = A.GetN();
596  int lwork = 2 * max(m, n);
597 
598 #ifdef SELDON_CHECK_DIMENSIONS
599  if (tau.GetM() < min(m, n))
600  throw WrongDim("GetQ_FromLQ", "tau not large enough");
601 #endif
602 
603  int k = min(m, n);
604  Vector<complex<double>, VectFull, Allocator1> work(lwork);
605  zungqr_(&n, &k, &k, A.GetDataVoid(), &n, tau.GetDataVoid(),
606  work.GetDataVoid(), &lwork, &info.GetInfoRef());
607 
608  if (m > n)
609  A.Resize(n, n);
610  }
611 
612 
613  // GETQ_FROMQR //
615 
616 
618  // GETLQ //
619 
620 
621  /*** ColMajor ***/
622 
623 
624  template<class Prop0, class Allocator0,
625  class Allocator1>
626  void GetLQ(Matrix<float, Prop0, ColMajor, Allocator0>& A,
627  Vector<float, VectFull, Allocator1>& tau,
628  LapackInfo& info)
629  {
630  int m = A.GetM();
631  int n = A.GetN();
632  int lwork = max(m,n);
633  Vector<float, VectFull, Allocator1> work(lwork);
634  tau.Reallocate(min(m, n));
635  sgelqf_(&m, &n, A.GetData(), &m, tau.GetData(),
636  work.GetData(), &lwork, &info.GetInfoRef());
637  }
638 
639 
640  template<class Prop0, class Allocator0,
641  class Allocator1>
642  void GetLQ(Matrix<double, Prop0, ColMajor, Allocator0>& A,
643  Vector<double, VectFull, Allocator1>& tau,
644  LapackInfo& info)
645  {
646  int m = A.GetM();
647  int n = A.GetN();
648  int lwork = max(m,n);
649  Vector<double, VectFull, Allocator1> work(lwork);
650  tau.Reallocate(min(m, n));
651  dgelqf_(&m, &n, A.GetData(), &m, tau.GetData(),
652  work.GetData(), &lwork, &info.GetInfoRef());
653  }
654 
655 
656  template<class Prop0, class Allocator0,
657  class Allocator1>
658  void GetLQ(Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
659  Vector<complex<float>, VectFull, Allocator1>& tau,
660  LapackInfo& info)
661  {
662  int m = A.GetM();
663  int n = A.GetN();
664  int lwork = max(m,n);
665  Vector<complex<float>, VectFull, Allocator1> work(lwork);
666  tau.Reallocate(min(m, n));
667  cgelqf_(&m, &n, A.GetData(), &m, tau.GetData(),
668  work.GetData(), &lwork, &info.GetInfoRef());
669  }
670 
671 
672  template<class Prop0, class Allocator0,
673  class Allocator1>
674  void GetLQ(Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
675  Vector<complex<double>, VectFull, Allocator1>& tau,
676  LapackInfo& info)
677  {
678  int m = A.GetM();
679  int n = A.GetN();
680  int lwork = max(m,n);
681  Vector<complex<double>, VectFull, Allocator1> work(lwork);
682  tau.Reallocate(min(m, n));
683  zgelqf_(&m, &n, A.GetData(), &m, tau.GetData(),
684  work.GetData(), &lwork, &info.GetInfoRef());
685  }
686 
687 
688  /*** RowMajor ***/
689 
690 
691  template<class Prop0, class Allocator0,
692  class Allocator1>
693  void GetLQ(Matrix<float, Prop0, RowMajor, Allocator0>& A,
694  Vector<float, VectFull, Allocator1>& tau,
695  LapackInfo& info)
696  {
697  int m = A.GetM();
698  int n = A.GetN();
699  int lwork = max(m,n);
700  Vector<float, VectFull, Allocator1> work(lwork);
701  tau.Reallocate(min(m, n));
702  // Factorization QR of A^t.
703  sgeqrf_(&n, &m, A.GetData(), &n, tau.GetData(),
704  work.GetData(), &lwork, &info.GetInfoRef());
705  }
706 
707 
708  template<class Prop0, class Allocator0,
709  class Allocator1>
710  void GetLQ(Matrix<double, Prop0, RowMajor, Allocator0>& A,
711  Vector<double, VectFull, Allocator1>& tau,
712  LapackInfo& info)
713  {
714  int m = A.GetM();
715  int n = A.GetN();
716  int lwork = max(m,n);
717  Vector<double, VectFull, Allocator1> work(lwork);
718  tau.Reallocate(min(m, n));
719  // Factorization LQ of A^t.
720  dgeqrf_(&n, &m, A.GetData(), &n, tau.GetData(),
721  work.GetData(), &lwork, &info.GetInfoRef());
722  }
723 
724 
725  template<class Prop0, class Allocator0,
726  class Allocator1>
727  void GetLQ(Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
728  Vector<complex<float>, VectFull, Allocator1>& tau,
729  LapackInfo& info)
730  {
731  int m = A.GetM();
732  int n = A.GetN();
733  int lwork = max(m,n);
734  Vector<complex<float>, VectFull, Allocator1> work(lwork);
735  tau.Reallocate(min(m, n));
736  // Factorization LQ of A^t.
737  cgeqrf_(&n, &m, A.GetData(), &n, tau.GetData(),
738  work.GetData(), &lwork, &info.GetInfoRef());
739  }
740 
741 
742  template<class Prop0, class Allocator0,
743  class Allocator1>
744  void GetLQ(Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
745  Vector<complex<double>, VectFull, Allocator1>& tau,
746  LapackInfo& info)
747  {
748  int m = A.GetM();
749  int n = A.GetN();
750  int lwork = max(m,n);
751  Vector<complex<double>, VectFull, Allocator1> work(lwork);
752  tau.Reallocate(min(m, n));
753  // Factorization LQ of A^t.
754  zgeqrf_(&n, &m, A.GetData(), &n, tau.GetData(),
755  work.GetData(), &lwork, &info.GetInfoRef());
756  }
757 
758 
759  // GETLQ //
761 
762 
764  // MLTQ_FROMQR //
765 
766 
767  /*** ColMajor ***/
768 
769 
770  template<class Prop0, class Allocator0,
771  class Allocator1, class Allocator2>
772  void MltQ_FromQR(const SeldonTranspose& trans,
773  const Matrix<float, Prop0, ColMajor, Allocator0>& A,
774  const Vector<float, VectFull, Allocator1>& tau,
775  Vector<float, VectFull, Allocator2>& b,
776  LapackInfo& info)
777  {
778  int lda = A.GetM();
779 
780 #ifdef SELDON_CHECK_DIMENSIONS
781  if (tau.GetM() < min(lda, A.GetN()))
782  throw WrongDim("MltQ_FromQR", "tau not large enough");
783 
784  if (b.GetM() < lda)
785  throw WrongDim("MltQ_FromQR", "b not large enough");
786 #endif
787 
788  int m = b.GetM();
789  int n = 1;
790  int k = tau.GetM();
791  int lwork = max(m, n);
792  Vector<float, VectFull, Allocator1> work(lwork);
793  char side('L');
794  char trans_ = trans.Char();
795  sormqr_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
796  b.GetData(), &m, work.GetData(), &lwork,
797  &info.GetInfoRef());
798  }
799 
800 
801  template<class Prop0, class Allocator0,
802  class Allocator1, class Allocator2>
803  void MltQ_FromQR(const SeldonTranspose& trans,
804  const Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
805  const Vector<complex<float>, VectFull, Allocator1>& tau,
806  Vector<complex<float>, VectFull, Allocator2>& b,
807  LapackInfo& info)
808  {
809  int lda = A.GetM();
810 
811 #ifdef SELDON_CHECK_DIMENSIONS
812  if (tau.GetM() < min(lda, A.GetN()))
813  throw WrongDim("MltQ_FromQR", "tau not large enough");
814 
815  if (b.GetM() < lda)
816  throw WrongDim("MltQ_FromQR", "b not large enough");
817 #endif
818 
819  int m = b.GetM();
820  int n = 1;
821  int k = tau.GetM();
822  int lwork = max(m, n);
823  Vector<complex<float>, VectFull, Allocator1> work(lwork);
824  char side('L');
825  char trans_ = trans.Char();
826  cunmqr_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
827  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
828  &info.GetInfoRef());
829  }
830 
831 
832  template<class Prop0, class Allocator0,
833  class Allocator1, class Allocator2>
834  void MltQ_FromQR(const SeldonTranspose& trans,
835  const Matrix<double, Prop0, ColMajor, Allocator0>& A,
836  const Vector<double, VectFull, Allocator1>& tau,
837  Vector<double, VectFull, Allocator2>& b,
838  LapackInfo& info)
839  {
840  int lda = A.GetM();
841 
842 #ifdef SELDON_CHECK_DIMENSIONS
843  if (tau.GetM() < min(lda, A.GetN()))
844  throw WrongDim("MltQ_FromQR", "tau not large enough");
845 
846  if (b.GetM() < lda)
847  throw WrongDim("MltQ_FromQR", "b not large enough");
848 #endif
849 
850  int m = b.GetM();
851  int n = 1;
852  int k = tau.GetM();
853  int lwork = max(m, n);
854  Vector<double, VectFull, Allocator1> work(lwork);
855  char side('L');
856  char trans_ = trans.Char();
857  dormqr_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
858  b.GetData(), &m, work.GetData(), &lwork,
859  &info.GetInfoRef());
860  }
861 
862 
863  template<class Prop0, class Allocator0,
864  class Allocator1, class Allocator2>
865  void MltQ_FromQR(const SeldonTranspose& trans,
866  const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
867  const Vector<complex<double>, VectFull, Allocator1>& tau,
868  Vector<complex<double>, VectFull, Allocator2>& b,
869  LapackInfo& info)
870  {
871  int lda = A.GetM();
872 
873 #ifdef SELDON_CHECK_DIMENSIONS
874  if (tau.GetM() < min(lda, A.GetN()))
875  throw WrongDim("MltQ_FromQR", "tau not large enough");
876 
877  if (b.GetM() < lda)
878  throw WrongDim("MltQ_FromQR", "b not large enough");
879 #endif
880 
881  int m = b.GetM();
882  int n = 1;
883  int k = tau.GetM();
884  int lwork = max(m, n);
885  Vector<complex<double>, VectFull, Allocator1> work(lwork);
886  char side('L');
887  char trans_ = trans.Char();
888  zunmqr_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
889  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
890  &info.GetInfoRef());
891  }
892 
893 
894  template<class Prop0, class Allocator0,
895  class Allocator1, class Allocator2>
896  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
897  const Matrix<float, Prop0, ColMajor, Allocator0>& A,
898  const Vector<float, VectFull, Allocator1>& tau,
899  Matrix<float, Prop0, ColMajor, Allocator2>& C,
900  LapackInfo& info)
901  {
902 
903 #ifdef SELDON_CHECK_DIMENSIONS
904  if (tau.GetM() < min(A.GetM(), A.GetN()))
905  throw WrongDim("MltQ_FromQR", "tau not large enough");
906 
907  if ( (side.Left() && C.GetM() < A.GetM())
908  || (side.Right() && C.GetN() < A.GetM()))
909  throw WrongDim("MltQ_FromQR", "C not large enough");
910 #endif
911 
912  int m = C.GetM();
913  int n = C.GetN();
914  int lwork = max(A.GetM(), A.GetN());
915  Vector<float, VectFull, Allocator1> work(lwork);
916  char side_ = side.Char();
917  char trans_ = trans.Char();
918  int lda = A.GetM();
919  int k = min(A.GetM(), A.GetN());
920 
921  sormqr_(&side_, &trans_, &m, &n, &k, A.GetData(), &lda,
922  tau.GetData(), C.GetData(), &m, work.GetData(), &lwork,
923  &info.GetInfoRef());
924  }
925 
926 
927  template<class Prop0, class Allocator0,
928  class Allocator1, class Allocator2>
929  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
930  const Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
931  const Vector<complex<float>, VectFull, Allocator1>& tau,
932  Matrix<complex<float>, Prop0, ColMajor, Allocator2>& C,
933  LapackInfo& info)
934  {
935 
936 #ifdef SELDON_CHECK_DIMENSIONS
937  if (tau.GetM() < min(A.GetM(), A.GetN()))
938  throw WrongDim("MltQ_FromQR", "tau not large enough");
939 
940  if ( (side.Left() && C.GetM() < A.GetM())
941  || (side.Right() && C.GetN() < A.GetM()))
942  throw WrongDim("MltQ_FromQR", "C not large enough");
943 #endif
944 
945  int m = C.GetM();
946  int n = C.GetN();
947  int lwork = max(A.GetM(), A.GetN());
948  Vector<complex<float>, VectFull, Allocator1> work(lwork);
949  char side_ = side.Char();
950  char trans_ = trans.Char();
951  int lda = A.GetM();
952  int k = min(A.GetM(), A.GetN());
953 
954  cunmqr_(&side_, &trans_, &m, &n, &k, A.GetDataVoid(), &lda,
955  tau.GetDataVoid(), C.GetDataVoid(), &m,
956  work.GetDataVoid(), &lwork, &info.GetInfoRef());
957  }
958 
959 
960  template<class Prop0, class Allocator0,
961  class Allocator1, class Allocator2>
962  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
963  const Matrix<double, Prop0, ColMajor, Allocator0>& A,
964  const Vector<double, VectFull, Allocator1>& tau,
965  Matrix<double, Prop0, ColMajor, Allocator2>& C,
966  LapackInfo& info)
967  {
968 
969 #ifdef SELDON_CHECK_DIMENSIONS
970  if (tau.GetM() < min(A.GetM(), A.GetN()))
971  throw WrongDim("MltQ_FromQR", "tau not large enough");
972 
973  if ( (side.Left() && C.GetM() < A.GetM())
974  || (side.Right() && C.GetN() < A.GetM()))
975  throw WrongDim("MltQ_FromQR", "C not large enough");
976 #endif
977 
978  int m = C.GetM();
979  int n = C.GetN();
980  int lwork = max(A.GetM(), A.GetN());
981  Vector<double, VectFull, Allocator1> work(lwork);
982  char side_ = side.Char();
983  char trans_ = trans.Char();
984  int lda = A.GetM();
985  int k = min(A.GetM(), A.GetN());
986 
987  dormqr_(&side_, &trans_, &m, &n, &k, A.GetData(), &lda,
988  tau.GetData(), C.GetData(), &m, work.GetData(), &lwork,
989  &info.GetInfoRef());
990  }
991 
992 
993  template<class Prop0, class Allocator0,
994  class Allocator1, class Allocator2>
995  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
996  const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
997  const Vector<complex<double>, VectFull, Allocator1>& tau,
998  Matrix<complex<double>, Prop0, ColMajor, Allocator2>& C,
999  LapackInfo& info)
1000  {
1001 
1002 #ifdef SELDON_CHECK_DIMENSIONS
1003  if (tau.GetM() < min(A.GetM(), A.GetN()))
1004  throw WrongDim("MltQ_FromQR", "tau not large enough");
1005 
1006  if ( (side.Left() && C.GetM() < A.GetM())
1007  || (side.Right() && C.GetN() < A.GetM()))
1008  throw WrongDim("MltQ_FromQR", "C not large enough");
1009 #endif
1010 
1011  int m = C.GetM();
1012  int n = C.GetN();
1013  int lwork = max(A.GetM(), A.GetN());
1014  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1015  char side_ = side.Char();
1016  char trans_ = trans.Char();
1017  int lda = A.GetM();
1018  int k = min(A.GetM(), A.GetN());
1019 
1020  zunmqr_(&side_, &trans_, &m, &n, &k, A.GetDataVoid(), &lda,
1021  tau.GetDataVoid(), C.GetDataVoid(), &m,
1022  work.GetDataVoid(), &lwork, &info.GetInfoRef());
1023  }
1024 
1025 
1026  template<class Prop0, class Allocator0,
1027  class Allocator1, class Allocator2>
1028  void MltQ_FromLQ(const SeldonTranspose& trans,
1029  const Matrix<float, Prop0, ColMajor, Allocator0>& A,
1030  const Vector<float, VectFull, Allocator1>& tau,
1031  Vector<float, VectFull, Allocator2>& b,
1032  LapackInfo& info)
1033  {
1034  int lda = A.GetM();
1035 
1036 #ifdef SELDON_CHECK_DIMENSIONS
1037  if (tau.GetM() < min(lda, A.GetN()))
1038  throw WrongDim("MltQ_FromQR", "tau not large enough");
1039 
1040  if (b.GetM() < A.GetN())
1041  throw WrongDim("MltQ_FromLQ", "b not large enough");
1042 #endif
1043 
1044  int m = b.GetM();
1045  int n = 1;
1046  int k = tau.GetM();
1047  int lwork = max(m, n);
1048  Vector<float, VectFull, Allocator1> work(lwork);
1049  char side('L');
1050  char trans_ = trans.Char();
1051  sormlq_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
1052  b.GetData(), &m, work.GetData(), &lwork,
1053  &info.GetInfoRef());
1054  }
1055 
1056 
1057  template<class Prop0, class Allocator0,
1058  class Allocator1, class Allocator2>
1059  void MltQ_FromLQ(const SeldonTranspose& trans,
1060  const Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
1061  const Vector<complex<float>, VectFull, Allocator1>& tau,
1062  Vector<complex<float>, VectFull, Allocator2>& b,
1063  LapackInfo& info)
1064  {
1065  int lda = A.GetM();
1066 
1067 #ifdef SELDON_CHECK_DIMENSIONS
1068  if (tau.GetM() < min(lda, A.GetN()))
1069  throw WrongDim("MltQ_FromQR", "tau not large enough");
1070 
1071  if (b.GetM() < A.GetN())
1072  throw WrongDim("MltQ_FromLQ", "b not large enough");
1073 #endif
1074 
1075  int m = b.GetM();
1076  int n = 1;
1077  int k = tau.GetM();
1078  int lwork = max(m, n);
1079  Vector<complex<float>, VectFull, Allocator1> work(lwork);
1080  char side('L');
1081  char trans_ = trans.Char();
1082  cunmlq_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
1083  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
1084  &info.GetInfoRef());
1085  }
1086 
1087 
1088  template<class Prop0, class Allocator0,
1089  class Allocator1, class Allocator2>
1090  void MltQ_FromLQ(const SeldonTranspose& trans,
1091  const Matrix<double, Prop0, ColMajor, Allocator0>& A,
1092  const Vector<double, VectFull, Allocator1>& tau,
1093  Vector<double, VectFull, Allocator2>& b,
1094  LapackInfo& info)
1095  {
1096  int lda = A.GetM();
1097 
1098 #ifdef SELDON_CHECK_DIMENSIONS
1099  if (tau.GetM() < min(lda, A.GetN()))
1100  throw WrongDim("MltQ_FromQR", "tau not large enough");
1101 
1102  if (b.GetM() < A.GetN())
1103  throw WrongDim("MltQ_FromLQ", "b not large enough");
1104 #endif
1105 
1106  int m = b.GetM();
1107  int n = 1;
1108  int k = tau.GetM();
1109  int lwork = max(m, n);
1110  Vector<double, VectFull, Allocator1> work(lwork);
1111  char side('L');
1112  char trans_ = trans.Char();
1113  dormlq_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
1114  b.GetData(), &m, work.GetData(), &lwork,
1115  &info.GetInfoRef());
1116  }
1117 
1118 
1119  template<class Prop0, class Allocator0,
1120  class Allocator1, class Allocator2>
1121  void MltQ_FromLQ(const SeldonTranspose& trans,
1122  const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
1123  const Vector<complex<double>, VectFull, Allocator1>& tau,
1124  Vector<complex<double>, VectFull, Allocator2>& b,
1125  LapackInfo& info)
1126  {
1127  int lda = A.GetM();
1128 
1129 #ifdef SELDON_CHECK_DIMENSIONS
1130  if (tau.GetM() < min(lda, A.GetN()))
1131  throw WrongDim("MltQ_FromQR", "tau not large enough");
1132 
1133  if (b.GetM() < A.GetN())
1134  throw WrongDim("MltQ_FromLQ", "b not large enough");
1135 #endif
1136 
1137  int m = b.GetM();
1138  int n = 1;
1139  int k = tau.GetM();
1140  int lwork = max(m, n);
1141  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1142  char side('L');
1143  char trans_ = trans.Char();
1144  zunmlq_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
1145  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
1146  &info.GetInfoRef());
1147  }
1148 
1149 
1150  template<class Prop0, class Allocator0,
1151  class Allocator1, class Allocator2>
1152  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1153  const Matrix<float, Prop0, ColMajor, Allocator0>& A,
1154  const Vector<float, VectFull, Allocator1>& tau,
1155  Matrix<float, Prop0, ColMajor, Allocator2>& C,
1156  LapackInfo& info)
1157  {
1158 
1159 #ifdef SELDON_CHECK_DIMENSIONS
1160  if (tau.GetM() < min(A.GetM(), A.GetN()))
1161  throw WrongDim("MltQ_FromQR", "tau not large enough");
1162 
1163  if ( (side.Left() && C.GetM() < A.GetN())
1164  || (side.Right() && C.GetN() < A.GetN()))
1165  throw WrongDim("MltQ_FromQR", "C not large enough");
1166 #endif
1167 
1168  int m = C.GetM();
1169  int n = C.GetN();
1170  int lwork = max(A.GetM(), A.GetN());
1171  Vector<float, VectFull, Allocator1> work(lwork);
1172  char side_ = side.Char();
1173  char trans_ = trans.Char();
1174  int lda = A.GetM();
1175  int k = min(A.GetM(), A.GetN());
1176 
1177  sormlq_(&side_, &trans_, &m, &n, &k, A.GetData(), &lda,
1178  tau.GetData(), C.GetData(), &m, work.GetData(), &lwork,
1179  &info.GetInfoRef());
1180  }
1181 
1182 
1183  template<class Prop0, class Allocator0,
1184  class Allocator1, class Allocator2>
1185  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1186  const Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
1187  const Vector<complex<float>, VectFull, Allocator1>& tau,
1188  Matrix<complex<float>, Prop0, ColMajor, Allocator2>& C,
1189  LapackInfo& info)
1190  {
1191 
1192 #ifdef SELDON_CHECK_DIMENSIONS
1193  if (tau.GetM() < min(A.GetM(), A.GetN()))
1194  throw WrongDim("MltQ_FromQR", "tau not large enough");
1195 
1196  if ( (side.Left() && C.GetM() < A.GetN())
1197  || (side.Right() && C.GetN() < A.GetN()))
1198  throw WrongDim("MltQ_FromQR", "C not large enough");
1199 #endif
1200 
1201  int m = C.GetM();
1202  int n = C.GetN();
1203  int lwork = max(A.GetM(), A.GetN());
1204  Vector<complex<float>, VectFull, Allocator1> work(lwork);
1205  char side_ = side.Char();
1206  char trans_ = trans.Char();
1207  int lda = A.GetM();
1208  int k = min(A.GetM(), A.GetN());
1209 
1210  cunmlq_(&side_, &trans_, &m, &n, &k, A.GetDataVoid(), &lda,
1211  tau.GetDataVoid(), C.GetDataVoid(), &m,
1212  work.GetDataVoid(), &lwork, &info.GetInfoRef());
1213  }
1214 
1215 
1216  template<class Prop0, class Allocator0,
1217  class Allocator1, class Allocator2>
1218  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1219  const Matrix<double, Prop0, ColMajor, Allocator0>& A,
1220  const Vector<double, VectFull, Allocator1>& tau,
1221  Matrix<double, Prop0, ColMajor, Allocator2>& C,
1222  LapackInfo& info)
1223  {
1224 
1225 #ifdef SELDON_CHECK_DIMENSIONS
1226  if (tau.GetM() < min(A.GetM(), A.GetN()))
1227  throw WrongDim("MltQ_FromQR", "tau not large enough");
1228 
1229  if ( (side.Left() && C.GetM() < A.GetN())
1230  || (side.Right() && C.GetN() < A.GetN()))
1231  throw WrongDim("MltQ_FromQR", "C not large enough");
1232 #endif
1233 
1234  int m = C.GetM();
1235  int n = C.GetN();
1236  int lwork = max(A.GetM(), A.GetN());
1237  Vector<double, VectFull, Allocator1> work(lwork);
1238  char side_ = side.Char();
1239  char trans_ = trans.Char();
1240  int lda = A.GetM();
1241  int k = min(A.GetM(), A.GetN());
1242 
1243  dormlq_(&side_, &trans_, &m, &n, &k, A.GetData(), &lda,
1244  tau.GetData(), C.GetData(), &m, work.GetData(), &lwork,
1245  &info.GetInfoRef());
1246  }
1247 
1248 
1249  template<class Prop0, class Allocator0,
1250  class Allocator1, class Allocator2>
1251  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1252  const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
1253  const Vector<complex<double>, VectFull, Allocator1>& tau,
1254  Matrix<complex<double>, Prop0, ColMajor, Allocator2>& C,
1255  LapackInfo& info)
1256  {
1257 
1258 #ifdef SELDON_CHECK_DIMENSIONS
1259  if (tau.GetM() < min(A.GetM(), A.GetN()))
1260  throw WrongDim("MltQ_FromQR", "tau not large enough");
1261 
1262  if ( (side.Left() && C.GetM() < A.GetN())
1263  || (side.Right() && C.GetN() < A.GetN()))
1264  throw WrongDim("MltQ_FromQR", "C not large enough");
1265 #endif
1266 
1267  int m = C.GetM();
1268  int n = C.GetN();
1269  int lwork = max(A.GetM(), A.GetN());
1270  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1271  char side_ = side.Char();
1272  char trans_ = trans.Char();
1273  int lda = A.GetM();
1274  int k = min(A.GetM(), A.GetN());
1275 
1276  zunmlq_(&side_, &trans_, &m, &n, &k, A.GetDataVoid(), &lda,
1277  tau.GetDataVoid(), C.GetDataVoid(), &m,
1278  work.GetDataVoid(), &lwork, &info.GetInfoRef());
1279  }
1280 
1281 
1282  /*** RowMajor ***/
1283 
1284 
1285  template<class Prop0, class Allocator0,
1286  class Allocator1, class Allocator2>
1287  void MltQ_FromQR(const SeldonTranspose& trans,
1288  const Matrix<float, Prop0, RowMajor, Allocator0>& A,
1289  const Vector<float, VectFull, Allocator1>& tau,
1290  Vector<float, VectFull, Allocator2>& b,
1291  LapackInfo& info)
1292  {
1293  int lda = A.GetN();
1294 
1295 #ifdef SELDON_CHECK_DIMENSIONS
1296  if (tau.GetM() < min(lda, A.GetM()))
1297  throw WrongDim("MltQ_FromQR", "tau not large enough");
1298 
1299  if (b.GetM() < A.GetM())
1300  throw WrongDim("MltQ_FromQR", "b not large enough");
1301 #endif
1302 
1303  int m = b.GetM();
1304  int n = 1;
1305  int k = tau.GetM();
1306  int lwork = max(m, n);
1307  Vector<float, VectFull, Allocator1> work(lwork);
1308  char side('L');
1309  char trans_ = trans.RevChar();
1310  sormlq_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
1311  b.GetData(), &m, work.GetData(), &lwork,
1312  &info.GetInfoRef());
1313  }
1314 
1315 
1316  template<class Prop0, class Allocator0,
1317  class Allocator1, class Allocator2>
1318  void MltQ_FromQR(const SeldonTranspose& trans,
1319  const Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
1320  const Vector<complex<float>, VectFull, Allocator1>& tau,
1321  Vector<complex<float>, VectFull, Allocator2>& b,
1322  LapackInfo& info)
1323  {
1324  int lda = A.GetN();
1325 
1326 #ifdef SELDON_CHECK_DIMENSIONS
1327  if (tau.GetM() < min(lda, A.GetM()))
1328  throw WrongDim("MltQ_FromQR", "tau not large enough");
1329 
1330  if (b.GetM() < A.GetM())
1331  throw WrongDim("MltQ_FromQR", "b not large enough");
1332 #endif
1333 
1334  int m = b.GetM();
1335  int n = 1;
1336  int k = tau.GetM();
1337  int lwork = max(m, n);
1338  Vector<complex<float>, VectFull, Allocator1> work(lwork);
1339  char side('L');
1340  char trans_('C');
1341  if (trans.ConjTrans())
1342  trans_ = 'N';
1343 
1344  Conjugate(b);
1345  cunmlq_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
1346  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
1347  &info.GetInfoRef());
1348 
1349  Conjugate(b);
1350  }
1351 
1352 
1353  template<class Prop0, class Allocator0,
1354  class Allocator1, class Allocator2>
1355  void MltQ_FromQR(const SeldonTranspose& trans,
1356  const Matrix<double, Prop0, RowMajor, Allocator0>& A,
1357  const Vector<double, VectFull, Allocator1>& tau,
1358  Vector<double, VectFull, Allocator2>& b,
1359  LapackInfo& info)
1360  {
1361  int lda = A.GetN();
1362 
1363 #ifdef SELDON_CHECK_DIMENSIONS
1364  if (tau.GetM() < min(lda, A.GetM()))
1365  throw WrongDim("MltQ_FromQR", "tau not large enough");
1366 
1367  if (b.GetM() < A.GetM())
1368  throw WrongDim("MltQ_FromQR", "b not large enough");
1369 #endif
1370 
1371  int m = b.GetM();
1372  int n = 1;
1373  int k = tau.GetM();
1374  int lwork = max(m, n);
1375  Vector<double, VectFull, Allocator1> work(lwork);
1376  char side('L');
1377  char trans_ = trans.RevChar();
1378  dormlq_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
1379  b.GetData(), &m, work.GetData(), &lwork,
1380  &info.GetInfoRef());
1381  }
1382 
1383 
1384  template<class Prop0, class Allocator0,
1385  class Allocator1, class Allocator2>
1386  void MltQ_FromQR(const SeldonTranspose& trans,
1387  const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
1388  const Vector<complex<double>, VectFull, Allocator1>& tau,
1389  Vector<complex<double>, VectFull, Allocator2>& b,
1390  LapackInfo& info)
1391  {
1392  int lda = A.GetN();
1393 
1394 #ifdef SELDON_CHECK_DIMENSIONS
1395  if (tau.GetM() < min(lda, A.GetM()))
1396  throw WrongDim("MltQ_FromQR", "tau not large enough");
1397 
1398  if (b.GetM() < A.GetM())
1399  throw WrongDim("MltQ_FromQR", "b not large enough");
1400 #endif
1401 
1402  int m = b.GetM();
1403  int n = 1;
1404  int k = tau.GetM();
1405  int lwork = max(m, n);
1406  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1407  char side('L');
1408  char trans_('C');
1409  if (trans.ConjTrans())
1410  trans_ = 'N';
1411 
1412  Conjugate(b);
1413  zunmlq_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
1414  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
1415  &info.GetInfoRef());
1416 
1417  Conjugate(b);
1418  }
1419 
1420 
1421  template<class Prop0, class Allocator0,
1422  class Allocator1, class Allocator2>
1423  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
1424  const Matrix<float, Prop0, RowMajor, Allocator0>& A,
1425  const Vector<float, VectFull, Allocator1>& tau,
1426  Matrix<float, Prop0, RowMajor, Allocator2>& C,
1427  LapackInfo& info)
1428  {
1429 
1430 #ifdef SELDON_CHECK_DIMENSIONS
1431  if (tau.GetM() < min(A.GetM(), A.GetN()))
1432  throw WrongDim("MltQ_FromQR", "tau not large enough");
1433 
1434  if ( (side.Left() && C.GetM() < A.GetM())
1435  || (side.Right() && C.GetN() < A.GetM()))
1436  throw WrongDim("MltQ_FromQR", "C not large enough");
1437 #endif
1438 
1439  int m = C.GetM();
1440  int n = C.GetN();
1441  int lwork = max(A.GetM(), A.GetN());
1442  Vector<float, VectFull, Allocator1> work(lwork);
1443  char side_ = side.RevChar();
1444  char trans_ = trans.Char();
1445  int lda = A.GetN();
1446  int k = min(A.GetM(), A.GetN());
1447 
1448  sormlq_(&side_, &trans_, &n, &m, &k, A.GetData(), &lda,
1449  tau.GetData(), C.GetData(), &n, work.GetData(), &lwork,
1450  &info.GetInfoRef());
1451  }
1452 
1453 
1454  template<class Prop0, class Allocator0,
1455  class Allocator1, class Allocator2>
1456  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
1457  const Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
1458  const Vector<complex<float>, VectFull, Allocator1>& tau,
1459  Matrix<complex<float>, Prop0, RowMajor, Allocator2>& C,
1460  LapackInfo& info)
1461  {
1462 
1463 #ifdef SELDON_CHECK_DIMENSIONS
1464  if (tau.GetM() < min(A.GetM(), A.GetN()))
1465  throw WrongDim("MltQ_FromQR", "tau not large enough");
1466 
1467  if ( (side.Left() && C.GetM() < A.GetM())
1468  || (side.Right() && C.GetN() < A.GetM()))
1469  throw WrongDim("MltQ_FromQR", "C not large enough");
1470 #endif
1471 
1472  int m = C.GetM();
1473  int n = C.GetN();
1474  int lwork = max(A.GetM(), A.GetN());
1475  Vector<complex<float>, VectFull, Allocator1> work(lwork);
1476  char side_ = side.RevChar();
1477  char trans_ = trans.Char();
1478  int lda = A.GetN();
1479  int k = min(A.GetM(), A.GetN());
1480 
1481  cunmlq_(&side_, &trans_, &n, &m, &k, A.GetDataVoid(), &lda,
1482  tau.GetDataVoid(), C.GetDataVoid(), &n,
1483  work.GetDataVoid(), &lwork, &info.GetInfoRef());
1484  }
1485 
1486 
1487  template<class Prop0, class Allocator0,
1488  class Allocator1, class Allocator2>
1489  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
1490  const Matrix<double, Prop0, RowMajor, Allocator0>& A,
1491  const Vector<double, VectFull, Allocator1>& tau,
1492  Matrix<double, Prop0, RowMajor, Allocator2>& C,
1493  LapackInfo& info)
1494  {
1495 
1496 #ifdef SELDON_CHECK_DIMENSIONS
1497  if (tau.GetM() < min(A.GetM(), A.GetN()))
1498  throw WrongDim("MltQ_FromQR", "tau not large enough");
1499 
1500  if ( (side.Left() && C.GetM() < A.GetM())
1501  || (side.Right() && C.GetN() < A.GetM()))
1502  throw WrongDim("MltQ_FromQR", "C not large enough");
1503 #endif
1504 
1505  int m = C.GetM();
1506  int n = C.GetN();
1507  int lwork = max(A.GetM(), A.GetN());
1508  Vector<double, VectFull, Allocator1> work(lwork);
1509  char side_ = side.RevChar();
1510  char trans_ = trans.Char();
1511  int lda = A.GetN();
1512  int k = min(A.GetM(), A.GetN());
1513 
1514  dormlq_(&side_, &trans_, &n, &m, &k, A.GetData(), &lda,
1515  tau.GetData(), C.GetData(), &n, work.GetData(), &lwork,
1516  &info.GetInfoRef());
1517  }
1518 
1519 
1520  template<class Prop0, class Allocator0,
1521  class Allocator1, class Allocator2>
1522  void MltQ_FromQR(const SeldonSide& side, const SeldonTranspose& trans,
1523  const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
1524  const Vector<complex<double>, VectFull, Allocator1>& tau,
1525  Matrix<complex<double>, Prop0, RowMajor, Allocator2>& C,
1526  LapackInfo& info)
1527  {
1528 
1529 #ifdef SELDON_CHECK_DIMENSIONS
1530  if (tau.GetM() < min(A.GetM(), A.GetN()))
1531  throw WrongDim("MltQ_FromQR", "tau not large enough");
1532 
1533  if ( (side.Left() && C.GetM() < A.GetM())
1534  || (side.Right() && C.GetN() < A.GetM()))
1535  throw WrongDim("MltQ_FromQR", "C not large enough");
1536 #endif
1537 
1538  int m = C.GetM();
1539  int n = C.GetN();
1540  int lwork = max(A.GetM(), A.GetN());
1541  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1542  char side_ = side.RevChar();
1543  char trans_ = trans.Char();
1544  int lda = A.GetN();
1545  int k = min(A.GetM(), A.GetN());
1546 
1547  zunmlq_(&side_, &trans_, &n, &m, &k, A.GetDataVoid(), &lda,
1548  tau.GetDataVoid(), C.GetDataVoid(), &n,
1549  work.GetDataVoid(), &lwork, &info.GetInfoRef());
1550  }
1551 
1552 
1553  template<class Prop0, class Allocator0,
1554  class Allocator1, class Allocator2>
1555  void MltQ_FromLQ(const SeldonTranspose& trans,
1556  const Matrix<float, Prop0, RowMajor, Allocator0>& A,
1557  const Vector<float, VectFull, Allocator1>& tau,
1558  Vector<float, VectFull, Allocator2>& b,
1559  LapackInfo& info)
1560  {
1561  int lda = A.GetN();
1562 
1563 #ifdef SELDON_CHECK_DIMENSIONS
1564  if (tau.GetM() < min(lda, A.GetM()))
1565  throw WrongDim("MltQ_FromLQ", "tau not large enough");
1566 
1567  if (b.GetM() < A.GetN())
1568  throw WrongDim("MltQ_FromLQ", "b not large enough");
1569 #endif
1570 
1571  int m = b.GetM();
1572  int n = 1;
1573  int k = tau.GetM();
1574  int lwork = max(m, n);
1575  Vector<float, VectFull, Allocator1> work(lwork);
1576  char side('L');
1577  char trans_ = trans.RevChar();
1578  sormqr_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
1579  b.GetData(), &m, work.GetData(), &lwork,
1580  &info.GetInfoRef());
1581  }
1582 
1583 
1584  template<class Prop0, class Allocator0,
1585  class Allocator1, class Allocator2>
1586  void MltQ_FromLQ(const SeldonTranspose& trans,
1587  const Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
1588  const Vector<complex<float>, VectFull, Allocator1>& tau,
1589  Vector<complex<float>, VectFull, Allocator2>& b,
1590  LapackInfo& info)
1591  {
1592  int lda = A.GetN();
1593 
1594 #ifdef SELDON_CHECK_DIMENSIONS
1595  if (tau.GetM() < min(lda, A.GetM()))
1596  throw WrongDim("MltQ_FromQR", "tau not large enough");
1597 
1598  if (b.GetM() < A.GetN())
1599  throw WrongDim("MltQ_FromLQ", "b not large enough");
1600 #endif
1601 
1602  int m = b.GetM();
1603  int n = 1;
1604  int k = tau.GetM();
1605  int lwork = max(m, n);
1606  Vector<complex<float>, VectFull, Allocator1> work(lwork);
1607  char side('L');
1608  char trans_('C');
1609  if (trans.ConjTrans())
1610  trans_ = 'N';
1611 
1612  Conjugate(b);
1613  cunmqr_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
1614  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
1615  &info.GetInfoRef());
1616 
1617  Conjugate(b);
1618  }
1619 
1620 
1621  template<class Prop0, class Allocator0,
1622  class Allocator1, class Allocator2>
1623  void MltQ_FromLQ(const SeldonTranspose& trans,
1624  const Matrix<double, Prop0, RowMajor, Allocator0>& A,
1625  const Vector<double, VectFull, Allocator1>& tau,
1626  Vector<double, VectFull, Allocator2>& b,
1627  LapackInfo& info)
1628  {
1629  int lda = A.GetN();
1630 
1631 #ifdef SELDON_CHECK_DIMENSIONS
1632  if (tau.GetM() < min(lda, A.GetM()))
1633  throw WrongDim("MltQ_FromLQ", "tau not large enough");
1634 
1635  if (b.GetM() < A.GetN())
1636  throw WrongDim("MltQ_FromLQ", "b not large enough");
1637 #endif
1638 
1639  int m = b.GetM();
1640  int n = 1;
1641  int k = tau.GetM();
1642  int lwork = max(m, n);
1643  Vector<double, VectFull, Allocator1> work(lwork);
1644  char side('L');
1645  char trans_ = trans.RevChar();
1646  dormqr_(&side, &trans_, &m, &n, &k, A.GetData(), &lda, tau.GetData(),
1647  b.GetData(), &m, work.GetData(), &lwork,
1648  &info.GetInfoRef());
1649  }
1650 
1651 
1652  template<class Prop0, class Allocator0,
1653  class Allocator1, class Allocator2>
1654  void MltQ_FromLQ(const SeldonTranspose& trans,
1655  const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
1656  const Vector<complex<double>, VectFull, Allocator1>& tau,
1657  Vector<complex<double>, VectFull, Allocator2>& b,
1658  LapackInfo& info)
1659  {
1660  int lda = A.GetN();
1661 
1662 #ifdef SELDON_CHECK_DIMENSIONS
1663  if (tau.GetM() < min(lda, A.GetM()))
1664  throw WrongDim("MltQ_FromQR", "tau not large enough");
1665 
1666  if (b.GetM() < A.GetN())
1667  throw WrongDim("MltQ_FromLQ", "b not large enough");
1668 #endif
1669 
1670  int m = b.GetM();
1671  int n = 1;
1672  int k = tau.GetM();
1673  int lwork = max(m, n);
1674  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1675  char side('L');
1676  char trans_('C');
1677  if (trans.ConjTrans())
1678  trans_ = 'N';
1679 
1680  Conjugate(b);
1681  zunmqr_(&side, &trans_, &m, &n, &k, A.GetDataVoid(), &lda, tau.GetDataVoid(),
1682  b.GetDataVoid(), &m, work.GetDataVoid(), &lwork,
1683  &info.GetInfoRef());
1684 
1685  Conjugate(b);
1686  }
1687 
1688 
1689  template<class Prop0, class Allocator0,
1690  class Allocator1, class Allocator2>
1691  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1692  const Matrix<float, Prop0, RowMajor, Allocator0>& A,
1693  const Vector<float, VectFull, Allocator1>& tau,
1694  Matrix<float, Prop0, RowMajor, Allocator2>& C,
1695  LapackInfo& info)
1696  {
1697 
1698 #ifdef SELDON_CHECK_DIMENSIONS
1699  if (tau.GetM() < min(A.GetM(), A.GetN()))
1700  throw WrongDim("MltQ_FromQR", "tau not large enough");
1701 
1702  if ( (side.Left() && C.GetM() < A.GetN())
1703  || (side.Right() && C.GetN() < A.GetN()))
1704  throw WrongDim("MltQ_FromQR", "C not large enough");
1705 #endif
1706 
1707  int m = C.GetM();
1708  int n = C.GetN();
1709  int lwork = max(A.GetM(), A.GetN());
1710  Vector<float, VectFull, Allocator1> work(lwork);
1711  char side_ = side.RevChar();
1712  char trans_ = trans.Char();
1713  int lda = A.GetN();
1714  int k = min(A.GetM(), A.GetN());
1715 
1716  sormqr_(&side_, &trans_, &n, &m, &k, A.GetData(), &lda,
1717  tau.GetData(), C.GetData(), &n, work.GetData(), &lwork,
1718  &info.GetInfoRef());
1719  }
1720 
1721 
1722  template<class Prop0, class Allocator0,
1723  class Allocator1, class Allocator2>
1724  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1725  const Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
1726  const Vector<complex<float>, VectFull, Allocator1>& tau,
1727  Matrix<complex<float>, Prop0, RowMajor, Allocator2>& C,
1728  LapackInfo& info)
1729  {
1730 #ifdef SELDON_CHECK_DIMENSIONS
1731  if (tau.GetM() < min(A.GetM(), A.GetN()))
1732  throw WrongDim("MltQ_FromQR", "tau not large enough");
1733 
1734  if ( (side.Left() && C.GetM() < A.GetN())
1735  || (side.Right() && C.GetN() < A.GetN()))
1736  throw WrongDim("MltQ_FromQR", "C not large enough");
1737 #endif
1738 
1739  int m = C.GetM();
1740  int n = C.GetN();
1741  int lwork = max(A.GetM(), A.GetN());
1742  Vector<complex<float>, VectFull, Allocator1> work(lwork);
1743  char side_ = side.RevChar();
1744  char trans_ = trans.Char();
1745  int lda = A.GetN();
1746  int k = min(A.GetM(), A.GetN());
1747 
1748  cunmqr_(&side_, &trans_, &n, &m, &k, A.GetDataVoid(), &lda,
1749  tau.GetDataVoid(), C.GetDataVoid(), &n,
1750  work.GetDataVoid(), &lwork, &info.GetInfoRef());
1751  }
1752 
1753 
1754  template<class Prop0, class Allocator0,
1755  class Allocator1, class Allocator2>
1756  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1757  const Matrix<double, Prop0, RowMajor, Allocator0>& A,
1758  const Vector<double, VectFull, Allocator1>& tau,
1759  Matrix<double, Prop0, RowMajor, Allocator2>& C,
1760  LapackInfo& info)
1761  {
1762 
1763 #ifdef SELDON_CHECK_DIMENSIONS
1764  if (tau.GetM() < min(A.GetM(), A.GetN()))
1765  throw WrongDim("MltQ_FromQR", "tau not large enough");
1766 
1767  if ( (side.Left() && C.GetM() < A.GetN())
1768  || (side.Right() && C.GetN() < A.GetN()))
1769  throw WrongDim("MltQ_FromQR", "C not large enough");
1770 #endif
1771 
1772  int m = C.GetM();
1773  int n = C.GetN();
1774  int lwork = max(A.GetM(), A.GetN());
1775  Vector<double, VectFull, Allocator1> work(lwork);
1776  char side_ = side.RevChar();
1777  char trans_ = trans.Char();
1778  int lda = A.GetN();
1779  int k = min(A.GetM(), A.GetN());
1780 
1781  dormqr_(&side_, &trans_, &n, &m, &k, A.GetData(), &lda,
1782  tau.GetData(), C.GetData(), &n, work.GetData(), &lwork,
1783  &info.GetInfoRef());
1784  }
1785 
1786 
1787  template<class Prop0, class Allocator0,
1788  class Allocator1, class Allocator2>
1789  void MltQ_FromLQ(const SeldonSide& side, const SeldonTranspose& trans,
1790  const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
1791  const Vector<complex<double>, VectFull, Allocator1>& tau,
1792  Matrix<complex<double>, Prop0, RowMajor, Allocator2>& C,
1793  LapackInfo& info)
1794  {
1795 #ifdef SELDON_CHECK_DIMENSIONS
1796  if (tau.GetM() < min(A.GetM(), A.GetN()))
1797  throw WrongDim("MltQ_FromQR", "tau not large enough");
1798 
1799  if ( (side.Left() && C.GetM() < A.GetN())
1800  || (side.Right() && C.GetN() < A.GetN()))
1801  throw WrongDim("MltQ_FromQR", "C not large enough");
1802 #endif
1803 
1804  int m = C.GetM();
1805  int n = C.GetN();
1806  int lwork = max(A.GetM(), A.GetN());
1807  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1808  char side_ = side.RevChar();
1809  char trans_ = trans.Char();
1810  int lda = A.GetN();
1811  int k = min(A.GetM(), A.GetN());
1812 
1813  zunmqr_(&side_, &trans_, &n, &m, &k, A.GetDataVoid(), &lda,
1814  tau.GetDataVoid(), C.GetDataVoid(), &n,
1815  work.GetDataVoid(), &lwork, &info.GetInfoRef());
1816  }
1817 
1818 
1819  // MLTQ_FROMQR //
1821 
1822 
1824  // SOLVEQR //
1825 
1826 
1827  /*** ColMajor ***/
1828 
1829 
1830  template <class Prop0, class Allocator0,
1831  class Allocator1,class Allocator2>
1832  void SolveQR(const Matrix<float, Prop0, ColMajor, Allocator0>& A,
1833  const Vector<float, VectFull, Allocator1>& tau,
1834  Vector<float, VectFull, Allocator2>& b,
1835  LapackInfo& info)
1836  {
1837  int m = A.GetM();
1838  int n = A.GetN();
1839  int k = tau.GetM();
1840 
1841  int nrhs = 1, nb = b.GetM();
1842  char side('L');
1843  char trans('T');
1844  int lwork = max(m, n);
1845  Vector<float, VectFull, Allocator1> work(lwork);
1846  // Computes Q^t b.
1847  sormqr_(&side, &trans, &m, &nrhs, &k, A.GetData(),
1848  &m, tau.GetData(), b.GetData(),
1849  &m, work.GetData(), &lwork, &info.GetInfoRef());
1850 
1851  b.Resize(n);
1852  for (int i = nb; i < n; i++)
1853  b(i) = 0;
1854 
1855  // Then solves R x = Q^t b.
1856  float alpha(1);
1857  cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans,
1858  CblasNonUnit, b.GetM(), nrhs,
1859  alpha, A.GetData(), A.GetM(), b.GetData(), b.GetM());
1860  }
1861 
1862 
1863  template <class Prop0, class Allocator0,
1864  class Allocator1,class Allocator2>
1865  void SolveQR(const Matrix<double, Prop0, ColMajor, Allocator0>& A,
1866  const Vector<double, VectFull, Allocator1>& tau,
1867  Vector<double, VectFull, Allocator2>& b,
1868  LapackInfo& info)
1869  {
1870  int m = A.GetM();
1871  int n = A.GetN();
1872  int k = tau.GetM();
1873 
1874 #ifdef SELDON_CHECK_DIMENSIONS
1875  if (tau.GetM() < min(m, n))
1876  throw WrongDim("SolveQR", "tau not large enough");
1877 
1878  if (b.GetM() < m)
1879  throw WrongDim("SolveQR", "b not large enough");
1880 #endif
1881 
1882  int nrhs = 1;
1883  char side('L');
1884  char trans('T');
1885  int lwork = max(m, n);
1886  Vector<double, VectFull, Allocator1> work(lwork);
1887  // Computes Q^t b.
1888  dormqr_(&side, &trans, &m, &nrhs, &k, A.GetData(),
1889  &m, tau.GetData(), b.GetData(),
1890  &lwork, work.GetData(), &lwork, &info.GetInfoRef());
1891 
1892  if (m > n)
1893  b.Resize(n);
1894 
1895  // Then solves R x = Q^t b.
1896  double alpha(1);
1897  cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans,
1898  CblasNonUnit, b.GetM(), nrhs,
1899  alpha, A.GetData(), A.GetM(), b.GetData(), b.GetM());
1900  }
1901 
1902 
1903  template <class Prop0, class Allocator0,
1904  class Allocator1,class Allocator2>
1905  void SolveQR(const Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
1906  const Vector<complex<float>, VectFull, Allocator1>& tau,
1907  Vector<complex<float>, VectFull, Allocator2>& b,
1908  LapackInfo& info)
1909  {
1910  int m = A.GetM();
1911  int n = A.GetN();
1912  int k = tau.GetM();
1913 
1914 #ifdef SELDON_CHECK_DIMENSIONS
1915  if (tau.GetM() < min(m, n))
1916  throw WrongDim("SolveQR", "tau not large enough");
1917 
1918  if (b.GetM() < m)
1919  throw WrongDim("SolveQR", "b not large enough");
1920 #endif
1921 
1922  int nrhs = 1;
1923  char side('L');
1924  char trans('C');
1925  int lwork = max(m, n);
1926  Vector<complex<float>, VectFull, Allocator1> work(lwork);
1927  // Computes Q^t b.
1928  cunmqr_(&side, &trans, &m, &nrhs, &k, A.GetData(),
1929  &m, tau.GetData(), b.GetData(),
1930  &m, work.GetData(), &lwork, &info.GetInfoRef());
1931 
1932  if (m > n)
1933  b.Resize(n);
1934 
1935  // Then solves R x = Q^t b.
1936  complex<float> alpha(1);
1937  cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans,
1938  CblasNonUnit, b.GetM(), nrhs,
1939  &alpha, A.GetData(), A.GetM(), b.GetData(), b.GetM());
1940  }
1941 
1942 
1943  template <class Prop0, class Allocator0,
1944  class Allocator1,class Allocator2>
1945  void SolveQR(const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
1946  const Vector<complex<double>, VectFull, Allocator1>& tau,
1947  Vector<complex<double>, VectFull, Allocator2>& b,
1948  LapackInfo& info)
1949  {
1950  int m = A.GetM();
1951  int n = A.GetN();
1952  int k = tau.GetM();
1953 
1954 #ifdef SELDON_CHECK_DIMENSIONS
1955  if (tau.GetM() < min(m, n))
1956  throw WrongDim("SolveQR", "tau not large enough");
1957 
1958  if (b.GetM() < m)
1959  throw WrongDim("SolveQR", "b not large enough");
1960 #endif
1961 
1962  int nrhs = 1;
1963  char side('L');
1964  char trans('C');
1965  int lwork = max(m, n);
1966  Vector<complex<double>, VectFull, Allocator1> work(lwork);
1967  // Computes Q^t b.
1968  zunmqr_(&side, &trans, &m, &nrhs, &k, A.GetData(),
1969  &m, tau.GetData(), b.GetData(),
1970  &m, work.GetData(), &lwork, &info.GetInfoRef());
1971 
1972  if (m > n)
1973  b.Resize(n);
1974 
1975  // Then solves R x = Q^t b.
1976  complex<double> alpha(1);
1977  cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans,
1978  CblasNonUnit, b.GetM(), nrhs,
1979  &alpha, A.GetData(), A.GetM(), b.GetData(), b.GetM());
1980  }
1981 
1982 
1983  /*** RowMajor ***/
1984 
1985 
1986  template <class Prop0, class Allocator0,
1987  class Allocator1,class Allocator2>
1988  void SolveQR(const Matrix<float, Prop0, RowMajor, Allocator0>& A,
1989  const Vector<float, VectFull, Allocator1>& tau,
1990  Vector<float, VectFull, Allocator2>& b,
1991  LapackInfo& info)
1992  {
1993  int m = A.GetM();
1994  int n = A.GetN();
1995 
1996 #ifdef SELDON_CHECK_DIMENSIONS
1997  if (tau.GetM() < min(m, n))
1998  throw WrongDim("SolveQR", "tau not large enough");
1999 
2000  if (b.GetM() < m)
2001  throw WrongDim("SolveQR", "b not large enough");
2002 #endif
2003 
2004  int k = tau.GetM();
2005  int nrhs = 1;
2006  char side('L');
2007  char trans('N');
2008  int lwork = max(m, n);
2009  Vector<float, VectFull, Allocator1> work(lwork);
2010  // Computes Q b.
2011  sormlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
2012  &n, tau.GetData(), b.GetData(),
2013  &m, work.GetData(), &lwork, &info.GetInfoRef());
2014 
2015  for (int i = m; i < n; i++)
2016  b(i) = 0;
2017 
2018  // Solves L^t y = b.
2019  float alpha(1);
2020  cblas_strsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
2021  CblasNonUnit, n, nrhs,
2022  alpha, A.GetData(), n, b.GetData(), b.GetM());
2023  }
2024 
2025 
2026  template <class Prop0, class Allocator0,
2027  class Allocator1,class Allocator2>
2028  void SolveQR(const Matrix<double, Prop0, RowMajor, Allocator0>& A,
2029  const Vector<double, VectFull, Allocator1>& tau,
2030  Vector<double, VectFull, Allocator2>& b,
2031  LapackInfo& info)
2032  {
2033  int m = A.GetM();
2034  int n = A.GetN();
2035 
2036 #ifdef SELDON_CHECK_DIMENSIONS
2037  if (tau.GetM() < min(m, n))
2038  throw WrongDim("SolveQR", "tau not large enough");
2039 
2040  if (b.GetM() < m)
2041  throw WrongDim("SolveQR", "b not large enough");
2042 #endif
2043 
2044  int k = min(m, n);
2045  int nrhs = 1;
2046  char side('L');
2047  char trans('N');
2048  int lwork = max(m, n);
2049  Vector<double, VectFull, Allocator1> work(lwork);
2050  // Computes Q b.
2051  dormlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
2052  &n, tau.GetData(), b.GetData(),
2053  &m, work.GetData(), &lwork, &info.GetInfoRef());
2054 
2055  if (m > n)
2056  b.Resize(n);
2057 
2058  // Solves L^t y = b.
2059  double alpha(1);
2060  cblas_dtrsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
2061  CblasNonUnit, k, nrhs,
2062  alpha, A.GetData(), n, b.GetData(), b.GetM());
2063 
2064 
2065  }
2066 
2067 
2068  template <class Prop0, class Allocator0,
2069  class Allocator1,class Allocator2>
2070  void SolveQR(const Matrix<double, Prop0, RowMajor, Allocator0>& A,
2071  const Vector<double, VectFull, Allocator1>& tau,
2072  Vector<complex<double>, VectFull, Allocator2>& b,
2073  LapackInfo& info)
2074  {
2075  int m = A.GetM();
2076  int n = A.GetN();
2077 
2078 #ifdef SELDON_CHECK_DIMENSIONS
2079  if (tau.GetM() < min(m, n))
2080  throw WrongDim("SolveQR", "tau not large enough");
2081 
2082  if (b.GetM() < m)
2083  throw WrongDim("SolveQR", "b not large enough");
2084 #endif
2085 
2086  int k = min(m, n);
2087  int nrhs = 2;
2088  char side('L');
2089  char trans('N');
2090  int lwork = max(m, n);
2091  Vector<double, VectFull, Allocator1> work(lwork);
2092  Vector<double> rhs(2*m);
2093  for (int i = 0; i < m; i++)
2094  {
2095  rhs(i) = realpart(b(i));
2096  rhs(i+m) = imagpart(b(i));
2097  }
2098 
2099  // Computes Q b
2100  dormlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
2101  &n, tau.GetData(), rhs.GetData(),
2102  &m, work.GetData(), &lwork, &info.GetInfoRef());
2103 
2104  for (int i = 0; i < n; i++)
2105  rhs(i+n) = rhs(i+m);
2106 
2107  rhs.Resize(2*n);
2108 
2109  // Solves L^t y = b.
2110  double alpha(1);
2111  cblas_dtrsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
2112  CblasNonUnit, k, nrhs,
2113  alpha, A.GetData(), n, rhs.GetData(), n);
2114 
2115  b.Reallocate(n);
2116  for (int i = 0; i < n; i++)
2117  b(i) = complex<double>(rhs(i), rhs(i+n));
2118  }
2119 
2120 
2121  template <class Prop0, class Allocator0,
2122  class Allocator1,class Allocator2>
2123  void SolveQR(const Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
2124  const Vector<complex<float>, VectFull, Allocator1>& tau,
2125  Vector<complex<float>, VectFull, Allocator2>& b,
2126  LapackInfo& info)
2127  {
2128  int m = A.GetM();
2129  int n = A.GetN();
2130  int k = tau.GetM();
2131 
2132 #ifdef SELDON_CHECK_DIMENSIONS
2133  if (tau.GetM() < min(m, n))
2134  throw WrongDim("SolveQR", "tau not large enough");
2135 
2136  if (b.GetM() < m)
2137  throw WrongDim("SolveQR", "b not large enough");
2138 #endif
2139 
2140  int nrhs = 1;
2141  char side('L');
2142  char trans('N');
2143  int lwork = max(m, n);
2144  Vector<complex<float>, VectFull, Allocator1> work(lwork);
2145  Conjugate(b);
2146  // Computes Q b.
2147  cunmlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
2148  &n, tau.GetData(), b.GetData(),
2149  &m, work.GetData(), &lwork, &info.GetInfoRef());
2150 
2151  Conjugate(b);
2152  if (m > n)
2153  b.Resize(n);
2154 
2155  // Solves L^t y = b.
2156  complex<float> alpha(1);
2157  cblas_ctrsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
2158  CblasNonUnit, k, nrhs,
2159  &alpha, A.GetData(), n, b.GetData(), b.GetM());
2160  }
2161 
2162 
2163  template <class Prop0, class Allocator0,
2164  class Allocator1,class Allocator2>
2165  void SolveQR(const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
2166  const Vector<complex<double>, VectFull, Allocator1>& tau,
2167  Vector<complex<double>, VectFull, Allocator2>& b,
2168  LapackInfo& info)
2169  {
2170  int m = A.GetM();
2171  int n = A.GetN();
2172  int k = tau.GetM();
2173 
2174 #ifdef SELDON_CHECK_DIMENSIONS
2175  if (tau.GetM() < min(m, n))
2176  throw WrongDim("SolveQR", "tau not large enough");
2177 
2178  if (b.GetM() < m)
2179  throw WrongDim("SolveQR", "b not large enough");
2180 #endif
2181 
2182  int nrhs = 1;
2183  char side('L');
2184  char trans('N');
2185  int lwork = max(m, n);
2186  Vector<complex<double>, VectFull, Allocator1> work(lwork);
2187  Conjugate(b);
2188  // Computes Q b.
2189  zunmlq_(&side, &trans, &m, &nrhs, &k, A.GetData(),
2190  &n, tau.GetData(), b.GetData(),
2191  &m, work.GetData(), &lwork, &info.GetInfoRef());
2192 
2193  Conjugate(b);
2194  if (m > n)
2195  b.Resize(n);
2196 
2197  // Solves L^t y = b.
2198  complex<double> alpha(1);
2199  cblas_ztrsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans,
2200  CblasNonUnit, k, nrhs,
2201  &alpha, A.GetData(), n, b.GetData(), b.GetM());
2202  }
2203 
2204 
2205  // SOLVEQR //
2207 
2208 
2210  // SOLVELQ //
2211 
2212 
2213  /*** ColMajor ***/
2214 
2215 
2216  template <class Prop0, class Allocator0,
2217  class Allocator1,class Allocator2>
2218  void SolveLQ(const Matrix<float, Prop0, ColMajor, Allocator0>& A,
2219  const Vector<float, VectFull, Allocator1>& tau,
2220  Vector<float, VectFull, Allocator2>& b,
2221  LapackInfo& info)
2222  {
2223  int m = A.GetM();
2224  int n = A.GetN();
2225  if (m > n)
2226  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2227 
2228  int k = tau.GetM();
2229  int nrhs = 1, nb = b.GetM();
2230  char side('L');
2231  char trans('T');
2232  int lwork = max(m, n);
2233  Vector<float, VectFull, Allocator1> work(lwork);
2234  // Solves L y = b.
2235  float alpha(1);
2236  cblas_strsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans,
2237  CblasNonUnit, m, nrhs,
2238  alpha, A.GetData(), m, b.GetData(), b.GetM());
2239 
2240  b.Resize(n);
2241  for (int i = nb; i < n; i++)
2242  b(i) = 0;
2243 
2244  // Computes Q^t b.
2245  sormlq_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2246  &m, tau.GetData(), b.GetData(),
2247  &n, work.GetData(), &lwork, &info.GetInfoRef());
2248  }
2249 
2250 
2251  template <class Prop0, class Allocator0,
2252  class Allocator1,class Allocator2>
2253  void SolveLQ(const Matrix<double, Prop0, ColMajor, Allocator0>& A,
2254  const Vector<double, VectFull, Allocator1>& tau,
2255  Vector<double, VectFull, Allocator2>& b,
2256  LapackInfo& info)
2257  {
2258  int m = A.GetM();
2259  int n = A.GetN();
2260  if (m > n)
2261  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2262 
2263 #ifdef SELDON_CHECK_DIMENSIONS
2264  if (tau.GetM() < min(m, n))
2265  throw WrongDim("SolveLQ", "tau not large enough");
2266 
2267  if (b.GetM() < n)
2268  throw WrongDim("SolveLQ", "b not large enough");
2269 #endif
2270 
2271  int k = tau.GetM();
2272  int nrhs = 1;
2273  char side('L');
2274  char trans('T');
2275  int lwork = max(m, n);
2276  Vector<double, VectFull, Allocator1> work(lwork);
2277  // Solves L y = b.
2278  double alpha(1);
2279  cblas_dtrsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans,
2280  CblasNonUnit, m, nrhs,
2281  alpha, A.GetData(), m, b.GetData(), b.GetM());
2282 
2283  // setting extra-components to 0
2284  for (int i = m; i < n; i++)
2285  b(i) = 0;
2286 
2287  // Applies Q^t b.
2288  dormlq_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2289  &m, tau.GetData(), b.GetData(),
2290  &n, work.GetData(), &lwork, &info.GetInfoRef());
2291  }
2292 
2293 
2294  template <class Prop0, class Allocator0,
2295  class Allocator1,class Allocator2>
2296  void SolveLQ(const Matrix<complex<float>, Prop0, ColMajor, Allocator0>& A,
2297  const Vector<complex<float>, VectFull, Allocator1>& tau,
2298  Vector<complex<float>, VectFull, Allocator2>& b,
2299  LapackInfo& info)
2300  {
2301  int m = A.GetM();
2302  int n = A.GetN();
2303  if (m > n)
2304  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2305 
2306 #ifdef SELDON_CHECK_DIMENSIONS
2307  if (tau.GetM() < min(m, n))
2308  throw WrongDim("SolveLQ", "tau not large enough");
2309 
2310  if (b.GetM() < n)
2311  throw WrongDim("SolveLQ", "b not large enough");
2312 #endif
2313 
2314  int k = tau.GetM();
2315  int nrhs = 1;
2316  char side('L');
2317  char trans('C');
2318  int lwork = max(m, n);
2319  Vector<complex<float>, VectFull, Allocator1> work(lwork);
2320  // Solve L y = b.
2321  complex<float> alpha(1);
2322  cblas_ctrsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans,
2323  CblasNonUnit, m, nrhs,
2324  &alpha, A.GetData(), m, b.GetData(), b.GetM());
2325 
2326  // setting extra-components to 0
2327  for (int i = m; i < n; i++)
2328  b(i) = 0;
2329 
2330  // Applies Q^t.
2331  cunmlq_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2332  &m, tau.GetData(), b.GetData(),
2333  &n, work.GetData(), &lwork, &info.GetInfoRef());
2334  }
2335 
2336 
2337  template <class Prop0, class Allocator0,
2338  class Allocator1,class Allocator2>
2339  void SolveLQ(const Matrix<complex<double>, Prop0, ColMajor, Allocator0>& A,
2340  const Vector<complex<double>, VectFull, Allocator1>& tau,
2341  Vector<complex<double>, VectFull, Allocator2>& b,
2342  LapackInfo& info)
2343  {
2344  int m = A.GetM();
2345  int n = A.GetN();
2346  if (m > n)
2347  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2348 
2349 #ifdef SELDON_CHECK_DIMENSIONS
2350  if (tau.GetM() < min(m, n))
2351  throw WrongDim("SolveLQ", "tau not large enough");
2352 
2353  if (b.GetM() < n)
2354  throw WrongDim("SolveLQ", "b not large enough");
2355 #endif
2356 
2357  int k = tau.GetM();
2358  int nrhs = 1;
2359  char side('L');
2360  char trans('C');
2361  int lwork = max(m, n);
2362  Vector<complex<double>, VectFull, Allocator1> work(lwork);
2363  // Solve L y = b.
2364  complex<double> alpha(1);
2365  cblas_ztrsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans,
2366  CblasNonUnit, m, nrhs,
2367  &alpha, A.GetData(), m, b.GetData(), b.GetM());
2368 
2369  // setting extra-components to 0
2370  for (int i = m; i < n; i++)
2371  b(i) = 0;
2372 
2373  // Applies Q^t.
2374  zunmlq_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2375  &m, tau.GetData(), b.GetData(),
2376  &n, work.GetData(), &lwork, &info.GetInfoRef());
2377  }
2378 
2379 
2380  /*** RowMajor ***/
2381 
2382 
2383  template <class Prop0, class Allocator0,
2384  class Allocator1,class Allocator2>
2385  void SolveLQ(const Matrix<float, Prop0, RowMajor, Allocator0>& A,
2386  const Vector<float, VectFull, Allocator1>& tau,
2387  Vector<float, VectFull, Allocator2>& b,
2388  LapackInfo& info)
2389  {
2390  int m = A.GetM();
2391  int n = A.GetN();
2392  if (m > n)
2393  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2394 
2395 #ifdef SELDON_CHECK_DIMENSIONS
2396  if (tau.GetM() < min(m, n))
2397  throw WrongDim("SolveLQ", "tau not large enough");
2398 
2399  if (b.GetM() < n)
2400  throw WrongDim("SolveLQ", "b not large enough");
2401 #endif
2402 
2403  int k = tau.GetM();
2404  int nrhs = 1;
2405  char side('L');
2406  char trans('N');
2407  int lwork = max(m, n);
2408  Vector<float, VectFull, Allocator1> work(lwork);
2409  // Solves R^t x = b.
2410  float alpha(1);
2411  cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,
2412  CblasNonUnit, k, nrhs,
2413  alpha, A.GetData(), A.GetN(), b.GetData(), b.GetM());
2414 
2415  for (int i = m; i < n; i++)
2416  b(i) = 0;
2417 
2418  // Multiplies by Q.
2419  sormqr_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2420  &n, tau.GetData(), b.GetData(),
2421  &n, work.GetData(), &lwork, &info.GetInfoRef());
2422  }
2423 
2424 
2425  template <class Prop0, class Allocator0,
2426  class Allocator1,class Allocator2>
2427  void SolveLQ(const Matrix<double, Prop0, RowMajor, Allocator0>& A,
2428  const Vector<double, VectFull, Allocator1>& tau,
2429  Vector<double, VectFull, Allocator2>& b,
2430  LapackInfo& info)
2431  {
2432  int m = A.GetM();
2433  int n = A.GetN();
2434  if (m > n)
2435  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2436 
2437 #ifdef SELDON_CHECK_DIMENSIONS
2438  if (tau.GetM() < min(m, n))
2439  throw WrongDim("SolveLQ", "tau not large enough");
2440 
2441  if (b.GetM() < n)
2442  throw WrongDim("SolveLQ", "b not large enough");
2443 #endif
2444 
2445  int k = tau.GetM();
2446  int nrhs = 1;
2447  char side('L');
2448  char trans('N');
2449  int lwork = max(m, n);
2450  Vector<double, VectFull, Allocator1> work(lwork);
2451  // Solves R^t x = b.
2452  double alpha(1);
2453  cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,
2454  CblasNonUnit, k, nrhs,
2455  alpha, A.GetData(), A.GetN(), b.GetData(), b.GetM());
2456 
2457  for (int i = m; i < n; i++)
2458  b(i) = 0;
2459 
2460  // Multiplies by Q.
2461  dormqr_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2462  &n, tau.GetData(), b.GetData(),
2463  &n, work.GetData(), &lwork, &info.GetInfoRef());
2464  }
2465 
2466 
2467  template <class Prop0, class Allocator0,
2468  class Allocator1,class Allocator2>
2469  void SolveLQ(const Matrix<complex<float>, Prop0, RowMajor, Allocator0>& A,
2470  const Vector<complex<float>, VectFull, Allocator1>& tau,
2471  Vector<complex<float>, VectFull, Allocator2>& b,
2472  LapackInfo& info)
2473  {
2474  int m = A.GetM();
2475  int n = A.GetN();
2476  if (m > n)
2477  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2478 
2479 #ifdef SELDON_CHECK_DIMENSIONS
2480  if (tau.GetM() < min(m, n))
2481  throw WrongDim("SolveLQ", "tau not large enough");
2482 
2483  if (b.GetM() < n)
2484  throw WrongDim("SolveLQ", "b not large enough");
2485 #endif
2486 
2487  int k = tau.GetM();
2488  int nrhs = 1;
2489  char side('L');
2490  char trans('N');
2491  int lwork = max(m, n);
2492  Vector<complex<float>, VectFull, Allocator1> work(lwork);
2493  // Solves R^t x = b.
2494  complex<float> alpha(1);
2495  cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,
2496  CblasNonUnit,k, nrhs,
2497  &alpha, A.GetData(), A.GetN(), b.GetData(), b.GetM());
2498 
2499  for (int i = m; i < n; i++)
2500  b(i) = 0;
2501 
2502  Conjugate(b);
2503  // Computes Q b.
2504  cunmqr_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2505  &n, tau.GetData(), b.GetData(),
2506  &n, work.GetData(), &lwork, &info.GetInfoRef());
2507 
2508  Conjugate(b);
2509  }
2510 
2511 
2512  template <class Prop0, class Allocator0,
2513  class Allocator1,class Allocator2>
2514  void SolveLQ(const Matrix<complex<double>, Prop0, RowMajor, Allocator0>& A,
2515  const Vector<complex<double>, VectFull, Allocator1>& tau,
2516  Vector<complex<double>, VectFull, Allocator2>& b,
2517  LapackInfo& info)
2518  {
2519  int m = A.GetM();
2520  int n = A.GetN();
2521  if (m > n)
2522  throw Undefined("SolveLQ", "Function implemented only for n >= m");
2523 
2524 #ifdef SELDON_CHECK_DIMENSIONS
2525  if (tau.GetM() < min(m, n))
2526  throw WrongDim("SolveLQ", "tau not large enough");
2527 
2528  if (b.GetM() < n)
2529  throw WrongDim("SolveLQ", "b not large enough");
2530 #endif
2531 
2532  int k = tau.GetM();
2533  int nrhs = 1;
2534  char side('L');
2535  char trans('N');
2536  int lwork = max(m, n);
2537  Vector<complex<double>, VectFull, Allocator1> work(lwork);
2538  // Solves R^t x = b.
2539  complex<double> alpha(1);
2540  cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans,
2541  CblasNonUnit,k, nrhs,
2542  &alpha, A.GetData(), A.GetN(), b.GetData(), b.GetM());
2543 
2544  for (int i = m; i < n; i++)
2545  b(i) = 0;
2546 
2547  Conjugate(b);
2548  // Computes Q b.
2549  zunmqr_(&side, &trans, &n, &nrhs, &k, A.GetData(),
2550  &n, tau.GetData(), b.GetData(),
2551  &n, work.GetData(), &lwork, &info.GetInfoRef());
2552 
2553  Conjugate(b);
2554  }
2555 
2556 
2557  // SOLVELQ //
2559 
2560 
2562 
2565  template<class T, class Prop, class Storage, class Allocator,
2566  class T2, class Prop2, class Storage2, class Allocator2>
2567  void GetLowerTriangular(const Matrix<T, Prop, Storage, Allocator>& A,
2568  Matrix<T2, Prop2, Storage2, Allocator2>& B)
2569  {
2570  int m = A.GetM();
2571  int n = A.GetN();
2572 
2573  T zero;
2574  SetComplexZero(zero);
2575 
2576  B.Clear();
2577  int k = min(m, n);
2578  B.Reallocate(m, k);
2579  B.Fill(zero);
2580 
2581  for (int i = 0; i < m; i++)
2582  for (int j = 0; j <= min(i, k-1); j++)
2583  B.Set(i, j, A(i, j));
2584 
2585  }
2586 
2587 
2589 
2592  template<class T, class Prop, class Storage, class Allocator,
2593  class T2, class Prop2, class Storage2, class Allocator2>
2594  void GetUpperTriangular(const Matrix<T, Prop, Storage, Allocator>& A,
2595  Matrix<T2, Prop2, Storage2, Allocator2>& B)
2596  {
2597  int m = A.GetM();
2598  int n = A.GetN();
2599 
2600  T zero;
2601  SetComplexZero(zero);
2602 
2603  B.Clear();
2604  B.Reallocate(min(m, n), n);
2605  B.Fill(zero);
2606 
2607  for (int i = 0; i < min(m, n); i++)
2608  for (int j = i; j < n; j++)
2609  B.Set(i, j, A(i, j));
2610 
2611  }
2612 
2613 
2615 
2618  template<class T, class Prop, class Storage, class Allocator,
2619  class T2, class Prop2, class Storage2, class Allocator2,
2620  class T3, class Prop3, class Storage3, class Allocator3>
2621  void GetQR(const Matrix<T, Prop, Storage, Allocator>& A,
2622  Matrix<T2, Prop2, Storage2, Allocator2>& Q,
2623  Matrix<T3, Prop3, Storage3, Allocator3>& R)
2624  {
2625  // QR-factorisation of A
2626  Vector<T> tau;
2627  Q = A;
2628  GetQR(Q, tau);
2629 
2630  // R is extracted
2631  GetUpperTriangular(Q, R);
2632 
2633  // then Q is explicited
2634  GetQ_FromQR(Q, tau);
2635  }
2636 
2637 
2639 
2642  template<class T, class Prop, class Storage, class Allocator,
2643  class T2, class Prop2, class Storage2, class Allocator2,
2644  class T3, class Prop3, class Storage3, class Allocator3>
2645  void GetLQ(const Matrix<T, Prop, Storage, Allocator>& A,
2646  Matrix<T2, Prop2, Storage2, Allocator2>& L,
2647  Matrix<T3, Prop3, Storage3, Allocator3>& Q)
2648  {
2649  // QR-factorisation of A
2650  Vector<T> tau;
2651  Q = A;
2652  GetLQ(Q, tau);
2653 
2654  // L is extracted
2655  GetLowerTriangular(Q, L);
2656 
2657  // then Q is explicited
2658  GetQ_FromLQ(Q, tau);
2659  }
2660 
2661 } // end namespace
2662 
2663 #define SELDON_FILE_LAPACK_LEAST_SQUARES_CXX
2664 #endif
Seldon::GetLowerTriangular
void GetLowerTriangular(const Matrix< T, Prop, Storage, Allocator > &A, Matrix< T2, Prop2, Storage2, Allocator2 > &B)
generic function to extract lower triangular part of a matrix
Definition: Lapack_LeastSquares.cxx:2567
Seldon::Vector< T >
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::Conjugate
void Conjugate(Matrix< T, Prop, Storage, Allocator > &A)
A is replaced by its conjugate.
Definition: Functions_Matrix.cxx:2915
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::GetUpperTriangular
void GetUpperTriangular(const Matrix< T, Prop, Storage, Allocator > &A, Matrix< T2, Prop2, Storage2, Allocator2 > &B)
generic function to extract upper triangular part of a matrix
Definition: Lapack_LeastSquares.cxx:2594