Home
Browse Files
  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        SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
 36       $                   INFO )
 37  *
 38  *  -- LAPACK routine (version 3.2) --
 39  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 40  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 41  *     November 2006
 42  */
 43 
 44 #ifndef FLENS_LAPACK_EIG_TREXC_TCC
 45 #define FLENS_LAPACK_EIG_TREXC_TCC 1
 46 
 47 #include <flens/blas/blas.h>
 48 #include <flens/lapack/lapack.h>
 49 
 50 namespace flens { namespace lapack {
 51 
 52 //== generic lapack implementation =============================================
 53 
 54 template <typename MT, typename MQ, typename IndexType, typename VWORK>
 55 IndexType
 56 trexc_generic(bool                          computeQ,
 57               GeMatrix<MT>                  &T,
 58               GeMatrix<MQ>                  &Q,
 59               IndexType                     &iFirst,
 60               IndexType                     &iLast,
 61               DenseVector<VWORK>            &work)
 62 {
 63     typedef typename GeMatrix<MT>::ElementType  ElementType;
 64 
 65     const IndexType     n = T.numRows();
 66     const ElementType   Zero(0);
 67 
 68     IndexType nBlockFirst, nBlockLast, nBlockNext, here;
 69     IndexType info = 0;
 70 //
 71 //  Quick return if possible
 72 //
 73     if (n<=1) {
 74         return info;
 75     }
 76 //
 77 //  Determine the first row of specified block
 78 //  and find out if it is 1 by 1 or 2 by 2.
 79 //
 80     if (iFirst>1) {
 81         if (T(iFirst,iFirst-1)!=Zero) {
 82             --iFirst;
 83         }
 84     }
 85     nBlockFirst = 1;
 86     if (iFirst<n) {
 87         if (T(iFirst+1,iFirst)!=Zero) {
 88             nBlockFirst = 2;
 89         }
 90     }
 91 //
 92 //  Determine the first row of the final block
 93 //  and find out if it is 1 by 1 or 2 by 2.
 94 //
 95     if (iLast>1) {
 96         if (T(iLast,iLast-1)!=Zero) {
 97             --iLast;
 98         }
 99     }
100     nBlockLast = 1;
101     if (iLast<n) {
102         if (T(iLast+1,iLast)!=Zero) {
103             nBlockLast = 2;
104         }
105     }
106 
107     if (iFirst==iLast) {
108         return info;
109     }
110 
111     if (iFirst<iLast) {
112 //
113 //      Update ILST
114 //
115         if (nBlockFirst==2 && nBlockLast==1) {
116             --iLast;
117         }
118         if (nBlockFirst==1 && nBlockLast==2) {
119             ++iLast;
120         }
121 
122         here = iFirst;
123 
124         do {
125 //
126 //          Swap block with next one below
127 //
128             if (nBlockFirst==1 || nBlockFirst==2) {
129 //
130 //              Current block either 1 by 1 or 2 by 2
131 //
132                 nBlockNext = 1;
133                 if (here+nBlockFirst+1<=n) {
134                     if (T(here+nBlockFirst+1,here+nBlockFirst)!=Zero) {
135                         nBlockNext = 2;
136                     }
137                 }
138                 info = laexc(computeQ, T, Q,
139                              here, nBlockFirst, nBlockNext,
140                              work);
141                 if (info!=0) {
142                     iLast = here;
143                     return info;
144                 }
145                 here += nBlockNext;
146 //
147 //              Test if 2 by 2 block breaks into two 1 by 1 blocks
148 //
149                 if (nBlockFirst==2) {
150                     if (T(here+1,here)==Zero) {
151                         nBlockFirst = 3;
152                     }
153                 }
154 
155             } else {
156 //
157 //              Current block consists of two 1 by 1 blocks each of which
158 //              must be swapped individually
159 //
160                 nBlockNext = 1;
161                 if (here+3<=n) {
162                     if (T(here+3,here+2)!=Zero) {
163                         nBlockNext = 2;
164                     }
165                 }
166                 info = laexc(computeQ, T, Q, here+1,
167                              IndexType(1), nBlockNext,
168                              work);
169                 if (info!=0) {
170                     iLast = here;
171                     return info;
172                 }
173                 if (nBlockNext==1) {
174 //
175 //                  Swap two 1 by 1 blocks, no problems possible
176 //
177                     laexc(computeQ, T, Q, here, IndexType(1), nBlockNext, work);
178                     ++here;
179                 } else {
180 //
181 //                  Recompute NBNEXT in case 2 by 2 split
182 //
183                     if (T(here+2,here+1)==Zero) {
184                         nBlockNext = 1;
185                     }
186                     if (nBlockNext==2) {
187 //
188 //                      2 by 2 Block did not split
189 //
190                         info = laexc(computeQ, T, Q, here,
191                                      IndexType(1), nBlockNext,
192                                      work);
193                         if (info!=0) {
194                             iLast = here;
195                             return info;
196                         }
197                         here += 2;
198                     } else {
199 //
200 //                      2 by 2 Block did split
201 //
202                         laexc(computeQ, T, Q,
203                               here, IndexType(1), IndexType(1),
204                               work);
205                         laexc(computeQ, T, Q,
206                               here+1, IndexType(1), IndexType(1),
207                               work);
208                         here += 2;
209                     }
210                 }
211             }
212         } while (here<iLast);
213     } else {
214 
215         here = iFirst;
216         do {
217 //
218 //          Swap block with next one above
219 //
220             if (nBlockFirst==1 || nBlockFirst==2) {
221 //
222 //              Current block either 1 by 1 or 2 by 2
223 //
224                 nBlockNext = 1;
225                 if (here>=3) {
226                     if (T(here-1,here-2)!=Zero) {
227                         nBlockNext = 2;
228                     }
229                 }
230                 info = laexc(computeQ, T, Q,
231                              here-nBlockNext, nBlockNext, nBlockFirst,
232                              work);
233                 if (info!=0) {
234                     iLast = here;
235                     return info;
236                 }
237                 here -= nBlockNext;
238 //
239 //              Test if 2 by 2 block breaks into two 1 by 1 blocks
240 //
241                 if (nBlockFirst==2) {
242                     if (T(here+1,here)==Zero) {
243                         nBlockFirst = 3;
244                     }
245                 }
246 
247             } else {
248 //
249 //              Current block consists of two 1 by 1 blocks each of which
250 //              must be swapped individually
251 //
252                 nBlockNext = 1;
253                 if (here>=3) {
254                     if (T(here-1,here-2)!=Zero) {
255                         nBlockNext = 2;
256                     }
257                 }
258                 info = laexc(computeQ, T, Q,
259                              here-nBlockNext, nBlockNext, IndexType(1),
260                              work);
261                 if (info!=0) {
262                     iLast = here;
263                     return info;
264                 } 
265                 if (nBlockNext==1) {
266 //
267 //                  Swap two 1 by 1 blocks, no problems possible
268 //
269                     laexc(computeQ, T, Q,
270                           here, nBlockNext, IndexType(1),
271                           work);
272                     --here;
273                 } else {
274 //
275 //                  Recompute NBNEXT in case 2 by 2 split
276 //
277                     if (T(here,here-1)==Zero) {
278                         nBlockNext = 1;
279                     }
280                     if (nBlockNext==2) {
281 //
282 //                      2 by 2 Block did not split
283 //
284                         info = laexc(computeQ, T, Q,
285                                      here-1, IndexType(2), IndexType(1),
286                                      work);
287                         if (info!=0) {
288                             iLast = here;
289                             return info;
290                         }
291                         here -= 2;
292                     } else {
293 //
294 //                      2 by 2 Block did split
295 //
296                         laexc(computeQ, T, Q,
297                               here, IndexType(1), IndexType(1),
298                               work);
299                         laexc(computeQ, T, Q,
300                               here-1, IndexType(1), IndexType(1),
301                               work);
302                         here -= 2;
303                     }
304                 }
305             }
306         } while (here>iLast);
307     }
308     iLast = here;
309     return info;
310 }
311 
312 //== interface for native lapack ===============================================
313 
314 #ifdef CHECK_CXXLAPACK
315 
316 template <typename MT, typename MQ, typename IndexType, typename VWORK>
317 IndexType
318 trexc_native(bool                          computeQ,
319              GeMatrix<MT>                  &T,
320              GeMatrix<MQ>                  &Q,
321              IndexType                     &iFirst,
322              IndexType                     &iLast,
323              DenseVector<VWORK>            &work)
324 {
325     typedef typename GeMatrix<MT>::ElementType ElementType;
326 
327     const char       COMPQ  = (computeQ) ? 'V' : 'N';
328     const INTEGER    N      = T.numRows();
329     const INTEGER    LDT    = T.leadingDimension();
330     const INTEGER    LDQ    = Q.leadingDimension();
331     INTEGER          IFST   = iFirst;
332     INTEGER          ILST   = iLast;
333     INTEGER          INFO;
334 
335     if (IsSame<ElementType,DOUBLE>::value) {
336         LAPACK_IMPL(dtrexc)(&COMPQ,
337                             &N,
338                             T.data(),
339                             &LDT,
340                             Q.data(),
341                             &LDQ,
342                             &IFST,
343                             &ILST,
344                             work.data(),
345                             &INFO);
346     } else {
347         ASSERT(0);
348     }
349     ASSERT(INFO>=0);
350 
351     iFirst  = IFST;
352     iLast   = ILST;
353     return INFO;
354 }
355 
356 #endif // CHECK_CXXLAPACK
357 
358 //== public interface ==========================================================
359 
360 template <typename MT, typename MQ, typename IndexType, typename VWORK>
361 IndexType
362 trexc(bool                          computeQ,
363       GeMatrix<MT>                  &T,
364       GeMatrix<MQ>                  &Q,
365       IndexType                     &iFirst,
366       IndexType                     &iLast,
367       DenseVector<VWORK>            &work)
368 {
369     LAPACK_DEBUG_OUT("trexc");
370 
371 //
372 //  Test the input parameters
373 //
374 #   ifndef NDEBUG
375     ASSERT(T.firstRow()==1);
376     ASSERT(T.firstCol()==1);
377     ASSERT(T.numRows()==T.numCols());
378 
379     const IndexType n = T.numRows();
380 
381     if (computeQ) {
382         ASSERT(Q.firstRow()==1);
383         ASSERT(Q.firstCol()==1);
384         ASSERT(Q.numRows()==Q.numCols());
385         ASSERT(Q.numRows()==n);
386     }
387 
388     ASSERT(work.firstIndex()==1);
389     ASSERT(work.length()==n);
390 #   endif
391 
392 #   ifdef CHECK_CXXLAPACK
393 //
394 //  Make copies of output arguments
395 //
396     typename GeMatrix<MT>::NoView           T_org      = T;
397     typename GeMatrix<MQ>::NoView           Q_org      = Q;
398     IndexType                               iFirst_org = iFirst;
399     IndexType                               iLast_org  = iLast;
400     typename DenseVector<VWORK>::NoView     work_org   = work;
401 #   endif
402 
403 //
404 //  Call implementation
405 //
406     IndexType info = trexc_generic(computeQ, T, Q, iFirst, iLast, work);
407 
408 #   ifdef CHECK_CXXLAPACK
409 //
410 //  Make copies of results computed by the generic implementation
411 //
412     typename GeMatrix<MT>::NoView           T_generic      = T;
413     typename GeMatrix<MQ>::NoView           Q_generic      = Q;
414     IndexType                               iFirst_generic = iFirst;
415     IndexType                               iLast_generic  = iLast;
416     typename DenseVector<VWORK>::NoView     work_generic   = work;
417 
418 //
419 //  restore output arguments
420 //
421     T       = T_org;
422     Q       = Q_org;
423     iFirst  = iFirst_org;
424     iLast   = iLast_org;
425     work    = work_org;
426 
427 //
428 //  Compare results
429 //
430     IndexType _info = trexc_native(computeQ, T, Q, iFirst, iLast, work);
431 
432     bool failed = false;
433     if (! isIdentical(T_generic, T, "T_generic", "T")) {
434         std::cerr << "CXXLAPACK: T_generic = " << T_generic << std::endl;
435         std::cerr << "F77LAPACK: T = " << T << std::endl;
436         failed = true;
437     }
438 
439     if (! isIdentical(Q_generic, Q, "Q_generic", "Q")) {
440         std::cerr << "CXXLAPACK: Q_generic = " << Q_generic << std::endl;
441         std::cerr << "F77LAPACK: Q = " << Q << std::endl;
442         failed = true;
443     }
444 
445     if (! isIdentical(iFirst_generic, iFirst, "iFirst_generic", "iFirst")) {
446         std::cerr << "CXXLAPACK:  iFirst_generic = "
447                   << iFirst_generic << std::endl;
448         std::cerr << "F77LAPACK: iFirst = "
449                   << iFirst << std::endl;
450         failed = true;
451     }
452 
453     if (! isIdentical(iLast_generic, iLast, "iLast_generic", "iLast")) {
454         std::cerr << "CXXLAPACK: iLast_generic = "
455                   << iLast_generic << std::endl;
456         std::cerr << "F77LAPACK: iLast = " << iLast << std::endl;
457         failed = true;
458     }
459 
460     if (! isIdentical(work_generic, work, "work_generic", "work")) {
461         std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
462         std::cerr << "F77LAPACK: work = " << work << std::endl;
463         failed = true;
464     }
465 
466     if (! isIdentical(info, _info, " info", "_info")) {
467         std::cerr << "CXXLAPACK:  info = " << info << std::endl;
468         std::cerr << "F77LAPACK: _info = " << _info << std::endl;
469         failed = true;
470     }
471 
472     if (failed) {
473         ASSERT(0);
474     }
475 #   endif
476 
477     return info;
478 }
479 
480 //-- forwarding ----------------------------------------------------------------
481 template <typename MT, typename MQ, typename IndexType, typename VWORK>
482 IndexType
483 trexc(bool                          computeQ,
484       MT                            &&T,
485       MQ                            &&Q,
486       IndexType                     &iFirst,
487       IndexType                     &iLast,
488       VWORK                         &&work)
489 {
490     CHECKPOINT_ENTER;
491     const IndexType info = trexc(computeQ, T, Q, iFirst, iLast, work);
492     CHECKPOINT_LEAVE;
493 
494     return info;
495 }
496 
497 } } // namespace lapack, flens
498 
499 #endif // FLENS_LAPACK_EIG_TREVC_TCC
Home