1 /*
  2  *   Copyright (c) 2011, Michael Lehn
  3  *
  4  *   All rights reserved.
  5  *
  6  *   Redistribution and use in source and binary forms, with or without
  7  *   modification, are permitted provided that the following conditions
  8  *   are met:
  9  *
 10  *   1) Redistributions of source code must retain the above copyright
 11  *      notice, this list of conditions and the following disclaimer.
 12  *   2) Redistributions in binary form must reproduce the above copyright
 13  *      notice, this list of conditions and the following disclaimer in
 14  *      the documentation and/or other materials provided with the
 15  *      distribution.
 16  *   3) Neither the name of the FLENS development group nor the names of
 17  *      its contributors may be used to endorse or promote products derived
 18  *      from this software without specific prior written permission.
 19  *
 20  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 21  *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 22  *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 23  *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 24  *   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 25  *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 26  *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 27  *   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 28  *   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 29  *   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 30  *   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 31  */
 32 
 33 /* Based on
 34  *
 35  *    INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 )
 36  *
 37  *    http://www.netlib.org/lapack/util/ilaenv.f
 38  *
 39  *  -- LAPACK auxiliary routine (version 3.2.1)                        --
 40  *
 41  *  -- April 2009                                                      --
 42  *
 43  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 44  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 45  *
 46  */
 47 
 48 #ifndef FLENS_LAPACK_AUX_ILAENV_TCC
 49 #define FLENS_LAPACK_AUX_ILAENV_TCC 1
 50 
 51 #include <complex>
 52 #include <string>
 53 #include <cstring>
 54 
 55 #include <flens/aux/issame.h>
 56 #include <flens/lapack/lapack.h>
 57 
 58 namespace flens { namespace lapack {
 59 
 60 //== generic lapack implementation =============================================
 61 
 62 template <int n>
 63 bool
 64 isSame(const char *a, const char *b)
 65 {
 66     for (int i=0; i<n; ++i) {
 67         if (a[i]!=b[i]) {
 68             return false;
 69         }
 70     }
 71     return true;
 72 }
 73 
 74 template <typename T>
 75 struct IsNotComplex
 76 {
 77     static const bool value = true;
 78 };
 79 
 80 template <typename T>
 81 struct IsNotComplex<std::complex<T> >
 82 {
 83     static const bool value = false;
 84 };
 85 
 86 template <typename T>
 87 struct IsComplex
 88 {
 89     static const bool value = !IsNotComplex<T>::value;
 90 };
 91 
 92 template <typename T>
 93 int
 94 ilaenv_generic(int spec, const char *name, const char *opts,
 95                int n1, int n2, int n3, int n4)
 96 {
 97     int result = -1;
 98     const char *c2 = name;
 99     const char *c3 = name + 2;
100     const char *c4 = name + 1;
101 
102     int nb, nbMin, nx;
103 
104     switch (spec) {
105 //
106 //      optimal blocksize
107 //
108         case 1:
109             nb = 32;
110             if (isSame<2>(c2, "GE")) {
111                 if (isSame<3>(c3, "TRF")) {
112                     if (IsNotComplex<T>::value) {
113                         nb = 64;
114                     } else {
115                         nb = 64;
116                     }
117                 } else if (isSame<3>(c3, "QRF") || isSame<3>(c3, "RQF")
118                        ||  isSame<3>(c3, "LQF") || isSame<3>(c3, "QLF")) {
119                     if (IsNotComplex<T>::value) {
120                         nb = 32;
121                     } else {
122                         nb = 32;
123                     }
124                 } else if (isSame<3>(c3, "HRD")) {
125                     if (IsNotComplex<T>::value) {
126                         nb = 32;
127                     } else {
128                         nb = 32;
129                     }
130                 } else if (isSame<3>(c3, "BRD")) {
131                     if (IsNotComplex<T>::value) {
132                         nb = 32;
133                     } else {
134                         nb = 32;
135                     }
136                 } else if (isSame<3>(c3, "TRI")) {
137                     if (IsNotComplex<T>::value) {
138                         nb = 64;
139                     } else {
140                         nb = 64;
141                     }
142                 }
143             } else if (isSame<2>(c2, "PO")) {
144                 if (isSame<3>(c3, "TRF")) {
145                     if (IsNotComplex<T>::value) {
146                         nb = 64;
147                     } else {
148                         nb = 64;
149                     }
150                 }
151             } else if (isSame<2>(c2, "SY")) {
152                 if (isSame<3>(c3, "TRF")) {
153                     if (IsNotComplex<T>::value) {
154                         nb = 64;
155                     } else {
156                         nb = 64;
157                     }
158                 } else if (IsNotComplex<T>::value && isSame<3>(c3, "TRD")) {
159                     nb = 32;
160                 } else if (IsNotComplex<T>::value && isSame<3>(c3, "GST")) {
161                     nb = 64;
162                 }
163             } else if (IsComplex<T>::value &&  isSame<2>(c2, "HE")) {
164                 if (isSame<3>(c3, "TRF")) {
165                     nb = 64;
166                 } else if (isSame<3>(c3, "TRD")) {
167                     nb = 32;
168                 } else if (isSame<3>(c3, "GST")) {
169                     nb = 64;
170                 }
171             } else if (IsNotComplex<T>::value && isSame<2>(c2, "OR")) {
172                 if (isSame<1>(c3, "G")) {
173                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
174                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
175                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
176                      || isSame<2>(c4, "BR"))
177                     {
178                         nb = 32;
179                     }
180                 } else if (isSame<1>(c3, "M")) {
181                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
182                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
183                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
184                      || isSame<2>(c4, "BR"))
185                     {
186                         nb = 32;
187                     }
188                 }
189             } else if (IsComplex<T>::value && isSame<2>(c2, "UN")) {
190                 if (isSame<1>(c3, "G")) {
191                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
192                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
193                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
194                      || isSame<2>(c4, "BR"))
195                     {
196                         nb = 32;
197                     }
198                } else if (isSame<1>(c3, "M")) {
199                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
200                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
201                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
202                      || isSame<2>(c4, "BR"))
203                     {
204                         nb = 32;
205                     }
206                }
207             } else if (isSame<2>(c2, "GB")) {
208                if (isSame<3>(c3, "TRF")) {
209                   if (IsNotComplex<T>::value) {
210                      if (n4<=64) {
211                         nb = 1;
212                      } else {
213                         nb = 32;
214                      }
215                   } else {
216                      if (n4<=64) {
217                         nb = 1;
218                      } else {
219                         nb = 32;
220                      }
221                   }
222                }
223             } else if (isSame<2>(c2, "PB")) {
224                if (isSame<3>(c3, "TRF")) {
225                   if (IsNotComplex<T>::value) {
226                      if (n2<=64) {
227                         nb = 1;
228                      } else {
229                         nb = 32;
230                      }
231                   } else {
232                      if (n2<=64 ) {
233                         nb = 1;
234                      } else {
235                         nb = 32;
236                      }
237                   }
238                }
239             } else if (isSame<2>(c2, "TR")) {
240                if (isSame<3>(c3, "TRI")) {
241                   if (IsNotComplex<T>::value) {
242                      nb = 64;
243                   } else {
244                      nb = 64;
245                   }
246                }
247             } else if (isSame<2>(c2, "LA")) {
248                if (isSame<3>(c3, "UUM")) {
249                   if (IsNotComplex<T>::value) {
250                      nb = 64;
251                   } else {
252                      nb = 64;
253                   }
254                }
255             } else if (IsNotComplex<T>::value && isSame<2>(c2, "ST") ) {
256                if (isSame<3>(c3, "EBZ")) {
257                   nb = 1;
258                }
259             }
260             result = nb*2;
261             break;
262 //
263 //      minimal blocksize
264 //
265         case 2:
266             nbMin = 2;
267             if (isSame<2>(c2, "GE")) {
268                 if (isSame<3>(c3, "QRF") || isSame<3>(c3, "RQF")
269                  || isSame<3>(c3, "LQF") || isSame<3>(c3, "QLF")) {
270                     if (IsNotComplex<T>::value ) {
271                          nbMin = 2;
272                     } else {
273                         nbMin = 2;
274                     }
275                 } else if (isSame<3>(c3, "HRD")) {
276                     if (IsNotComplex<T>::value ) {
277                         nbMin = 2;
278                     } else {
279                         nbMin = 2;
280                     }
281                 } else if (isSame<3>(c3, "BRD")) {
282                     if (IsNotComplex<T>::value ) {
283                         nbMin = 2;
284                     } else {
285                         nbMin = 2;
286                     }
287                 } else if (isSame<3>(c3, "TRI")) {
288                     if (IsNotComplex<T>::value ) {
289                         nbMin = 2;
290                     } else {
291                         nbMin = 2;
292                     }
293                 }
294             } else if (isSame<2>(c2, "SY")) {
295                 if (isSame<3>(c3, "TRF")) {
296                     if (IsNotComplex<T>::value ) {
297                         nbMin = 8;
298                     } else {
299                         nbMin = 8;
300                     }
301                 } else if (IsNotComplex<T>::value && isSame<3>(c3, "TRD")) {
302                     nbMin = 2;
303                 }
304             } else if (IsComplex<T>::value && isSame<2>(c2, "HE")) {
305                 if (isSame<3>(c3, "TRD")) {
306                     nbMin = 2;
307                 }
308             } else if (IsNotComplex<T>::value && isSame<2>(c2, "OR")) {
309                 if (isSame<1>(c3, "G")) {
310                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
311                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
312                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
313                      || isSame<2>(c4, "BR")) {
314                         nbMin = 2;
315                     }
316                 } else if (isSame<1>(c3, "M")) {
317                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
318                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
319                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
320                      || isSame<2>(c4, "BR")) {
321                         nbMin = 2;
322                     }
323                 }
324             } else if (IsComplex<T>::value && isSame<2>(c2, "UN")) {
325                 if (isSame<1>(c3, "G")) {
326                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
327                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
328                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
329                      || isSame<2>(c4, "BR")) {
330                         nbMin = 2;
331                     }
332                 } else if (isSame<1>(c3, "M")) {
333                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
334                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
335                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
336                      || isSame<2>(c4, "BR")) {
337                         nbMin = 2;
338                     }
339                 }
340             }
341             result = nbMin;
342             break;
343 //
344 //      crossover point
345 //
346         case 3:
347             nx = 32;
348             if (isSame<2>(c2, "GE")) {
349                 if (isSame<3>(c3, "QRF") || isSame<3>(c3, "RQF")
350                  || isSame<3>(c3, "LQF") || isSame<3>(c3, "QLF")) {
351                     if (IsNotComplex<T>::value) {
352                         nx =  128;
353                     } else {
354                         nx =  128;
355                     }
356                 } else if (isSame<3>(c3, "HRD")) {
357                     if (IsNotComplex<T>::value) {
358                         nx =  128;
359                     } else {
360                         nx =  128;
361                     }
362                 } else if (isSame<3>(c3, "BRD")) {
363                     if (IsNotComplex<T>::value) {
364                         nx =  128;
365                     } else {
366                         nx =  128;
367                     }
368                 }
369             } else if (isSame<2>(c2, "SY")) {
370                 if (IsNotComplex<T>::value && isSame<3>(c3, "TRD")) {
371                     nx =  32;
372                 }
373             } else if (IsComplex<T>::value && isSame<2>(c2, "HE")) {
374                 if (isSame<3>(c3, "TRD")) {
375                     nx =  32;
376                 }
377             } else if (IsNotComplex<T>::value && isSame<2>(c2, "OR")) {
378                 if (isSame<1>(c3, "G")) {
379                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
380                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
381                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
382                      || isSame<2>(c4, "BR")) {
383                         nx =  128;
384                     }
385                 }
386             } else if (IsComplex<T>::value && isSame<2>(c2, "UN")) {
387                 if (isSame<1>(c3, "G")) {
388                     if (isSame<2>(c4, "QR") || isSame<2>(c4, "RQ")
389                      || isSame<2>(c4, "LQ") || isSame<2>(c4, "QL")
390                      || isSame<2>(c4, "HR") || isSame<2>(c4, "TR")
391                      || isSame<2>(c4, "BR")) {
392                         nx =  128;
393                     }
394                 }
395             }
396             result = nx/1.5;
397             break;
398 //
399 //      DEPRECATED: number of shifts used in the nonsymmetric eigenvalue
400 //                  routines
401 //
402         case 4:
403             ASSERT(0);
404             break;
405 //
406 //      12 <= pec<= 16: hseqr or one of its subroutines .. 
407 //
408         case 12:
409         case 13:
410         case 14:
411         case 15:
412         case 16:
413             result = iparmq<T>(spec, name, opts, n1, n2, n3, n4);
414             break;
415 
416         default:
417             ASSERT(0);
418             result = -1;
419             break;
420 
421     }
422 
423 #   ifdef LOG_ILAENV
424     std::cerr << "ILAENV_GENERIC: "
425               << "spec = " << spec
426               << ", name = " << name
427               << ", opts = " << opts
428               << ", n1 " << n1
429               << ", n2 " << n2
430               << ", n3 " << n3
431               << ", n4 " << n4
432               << ", result = " << result
433               << std::endl;
434 #   endif
435 
436 
437     return result;
438 }
439 
440 //== interface for native lapack ===============================================
441 
442 #if defined CHECK_CXXLAPACK || defined USE_NATIVE_ILAENV
443 
444 #ifndef LAPACK_DECL
445 #   define  LAPACK_DECL(x)    x##_
446 
447 #   define  INTEGER           int
448 #endif
449 
450 extern "C" {
451 
452 INTEGER
453 LAPACK_DECL(ilaenv)(const INTEGER   *SPEC,
454                     const char      *NAME,
455                     const char      *OPTS,
456                     const INTEGER   *N1,
457                     const INTEGER   *N2,
458                     const INTEGER   *N3,
459                     const INTEGER   *N4,
460                     int             NAME_LEN,
461                     int             OPTS_LEN);
462 
463 // extern "C"
464 
465 
466 
467 template <typename T>
468 int
469 ilaenv_LapackTest(int spec, const char *_name, const char *_opts,
470                   int n1, int n2, int n3, int n4)
471 {
472     using std::string;
473     using std::complex;
474 
475     string opts(_opts);
476     string name;
477     if (IsSame<T,float>::value) {
478         name = string("S") + string(_name);
479     } else if (IsSame<T,double>::value) {
480         name = string("D") + string(_name);
481     } else if (IsSame<T,complex<float> >::value) {
482         name = string("C") + string(_name);
483     } else if (IsSame<T,complex<double> >::value) {
484         name = string("Z") + string(_name);
485     } else {
486         ASSERT(0);
487     }
488 
489 #if defined CHECK_CXXLAPACK || defined USE_NATIVE_ILAENV
490     int result = LAPACK_DECL(ilaenv)(&spec, name.c_str(), opts.c_str(),
491                                      &n1, &n2, &n3, &n4,
492                                      strlen(name.c_str()),
493                                      strlen(opts.c_str()));
494 
495 #   ifdef LOG_ILAENV
496     std::cerr << "ILAENV_GENERIC: "
497               << "spec = " << spec
498               << ", name = " << name
499               << ", opts = " << opts
500               << ", n1 " << n1
501               << ", n2 " << n2
502               << ", n3 " << n3
503               << ", n4 " << n4
504               << ", result = " << result
505               << std::endl;
506 #   endif
507 
508 
509     return result;
510 #else
511     ASSERT(0);
512 #endif
513 }
514 
515 #endif // CHECK_CXXLAPACK
516 
517 //== public interface ==========================================================
518 
519 template <typename T>
520 int
521 ilaenv(int spec, const char *name, const char *opts,
522        int n1, int n2, int n3, int n4)
523 {
524 #if defined CHECK_CXXLAPACK || defined USE_NATIVE_ILAENV
525     return ilaenv_LapackTest<T>(spec, name, opts, n1, n2, n3, n4);
526 #else
527     return ilaenv_generic<T>(spec, name, opts, n1, n2, n3, n4);
528 #endif
529 }
530 
531 } } // namespace lapack, flens
532 
533 #endif // FLENS_LAPACK_AUX_ILAENV_TCC