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 DLACN2( N, V, X, ISGN, EST, KASE, ISAVE )
 36  *
 37  *  -- LAPACK auxiliary routine (version 3.2) --
 38  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 39  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 40  *     November 2006
 41  *
 42  */
 43 
 44 #ifndef FLENS_LAPACK_AUX_LACN2_TCC
 45 #define FLENS_LAPACK_AUX_LACN2_TCC 1
 46 
 47 #include <cmath>
 48 #include <limits>
 49 
 50 namespace flens { namespace lapack {
 51 
 52 //== generic lapack implementation =============================================
 53 template <typename  VV, typename VX, typename VSGN, typename EST,
 54           typename IndexType, typename VSAVE>
 55 void
 56 lacn2_generic(DenseVector<VV> &v, DenseVector<VX> &x, DenseVector<VSGN> &sgn,
 57               EST &est, IndexType &kase, DenseVector<VSAVE> &iSave)
 58 {
 59     using std::abs;
 60 
 61     typedef typename DenseVector<VV>::ElementType    T;
 62     const T Zero(0), One(1), Two(2);
 63 
 64     const IndexType itMax = 5;
 65     const IndexType n = v.length();
 66 //  ..
 67 //  .. Executable Statements ..
 68 //
 69     if (kase==0) {
 70         for (IndexType i=1; i<=n; ++i) {
 71             x(i) = One / T(n);
 72         }
 73         kase = 1;
 74         iSave(1) = 1;
 75         return;
 76     }
 77 
 78     IndexType jLast;
 79     T         altSign, estOld;
 80     bool      converged;
 81 
 82     switch (iSave(1)) {
 83 //
 84 //  ................ ENTRY   (ISAVE( 1 ) = 1)
 85 //  FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
 86 //
 87     case 1:
 88         if (n==1) {
 89             v(1) = x(1);
 90             est = abs(v(1));
 91     //      ... QUIT
 92             goto QUIT;
 93         }
 94         est = blas::asum(x);
 95 
 96         for (IndexType i=1; i<=n; ++i) {
 97             x(i)    = sign(One, x(i));
 98             sgn(i)  = nint(x(i));
 99         }
100         kase = 2;
101         iSave(1) = 2;
102         return;
103 //
104 //  ................ ENTRY   (ISAVE( 1 ) = 2)
105 //  FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
106 //
107     case 2:
108         iSave(2) = blas::iamax(x);
109         iSave(3) = 2;
110 //
111 //      MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
112 //
113     MAIN:
114         for (IndexType i=1; i<=n; ++i) {
115             x(i) = Zero;
116         }
117         x(iSave(2)) = One;
118         kase = 1;
119         iSave(1) = 3;
120         return;
121 //
122 //  ................ ENTRY   (ISAVE( 1 ) = 3)
123 //  X HAS BEEN OVERWRITTEN BY A*X.
124 //
125     case 3:
126         v = x;
127         estOld = est;
128         est = blas::asum(v);
129 
130         converged = true;
131         for (IndexType i=1; i<=n; ++i) {
132             if (nint(sign(One,x(i)))!=sgn(i)) {
133                 converged = false;
134                 break;
135             }
136         }
137         if (converged) {
138 //          REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
139             goto ITERATION_COMPLETE;
140         }
141 
142 //      TEST FOR CYCLING.
143         if (est<=estOld) {
144             goto ITERATION_COMPLETE;
145         }
146 
147         for (IndexType i=1; i<=n; ++i) {
148             x(i)    = sign(One, x(i));
149             sgn(i)  = nint(x(i));
150         }
151         kase = 2;
152         iSave(1) = 4;
153         return;
154 //
155 //  ................ ENTRY   (ISAVE( 1 ) = 4)
156 //  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
157 //
158     case 4:
159         jLast = iSave(2);
160         iSave(2) = blas::iamax(x);
161         if (x(jLast)!=abs(x(iSave(2))) && iSave(3)<itMax) {
162             ++iSave(3);
163             goto MAIN;
164         }
165 //
166 //      ITERATION COMPLETE.  FINAL STAGE.
167 //
168     ITERATION_COMPLETE:
169         altSign = One;
170         for (IndexType i=1; i<=n; ++i) {
171             x(i) = altSign*(One + T(i-1)/T(n-1));
172             altSign = -altSign;
173         }
174         kase = 1;
175         iSave(1) = 5;
176         return;
177 //
178 //  ................ ENTRY   (ISAVE( 1 ) = 5)
179 //  X HAS BEEN OVERWRITTEN BY A*X.
180 //
181     case 5:
182         const T temp = Two*(blas::asum(x) / T(3*n));
183         if (temp>est) {
184             v = x;
185             est = temp;
186         }
187     }
188 
189 QUIT:
190     kase = 0;
191 }
192 
193 //== interface for native lapack ===============================================
194 
195 #ifdef CHECK_CXXLAPACK
196 
197 
198 template <typename  VV, typename VX, typename VSGN, typename EST,
199           typename IndexType, typename VSAVE>
200 void
201 lacn2_native(DenseVector<VV> &v, DenseVector<VX> &x, DenseVector<VSGN> &sgn,
202              EST &est, IndexType &kase, DenseVector<VSAVE> &iSave)
203 {
204     typedef typename DenseVector<VV>::ElementType   T;
205 
206     const INTEGER                   N = v.length();
207     DenseVector<Array<INTEGER> >    _sgn(N);
208     T                               _EST = est;
209     INTEGER                         _KASE = kase;
210     DenseVector<Array<INTEGER> >    _iSave(3);
211 
212     for (INTEGER i=1; i<=N; ++i) {
213         _sgn(i) = sgn(i);
214     }
215     for (INTEGER i=1; i<=3; ++i) {
216         _iSave(i) = iSave(i);
217     }
218     _iSave = iSave;
219 
220     if (IsSame<T,DOUBLE>::value) {
221         LAPACK_IMPL(dlacn2)(&N,
222                             v.data(),
223                             x.data(),
224                             _sgn.data(),
225                             &_EST,
226                             &_KASE,
227                             _iSave.data());
228     } else {
229         ASSERT(0);
230     }
231 
232     est  = _EST;
233     kase = _KASE;
234 
235     for (INTEGER i=1; i<=N; ++i) {
236         sgn(i) = _sgn(i);
237     }
238     for (INTEGER i=1; i<=3; ++i) {
239         iSave(i) = _iSave(i);
240     }
241 }
242 
243 #endif // CHECK_CXXLAPACK
244 
245 //== public interface ==========================================================
246 template <typename  VV, typename VX, typename VSGN, typename EST,
247           typename IndexType, typename VSAVE>
248 void
249 lacn2(DenseVector<VV> &v, DenseVector<VX> &x, DenseVector<VSGN> &sgn,
250       EST &est, IndexType &kase, DenseVector<VSAVE> &iSave)
251 {
252 //
253 //  Test the input parameters
254 //
255 #   ifndef NDEBUG
256     const IndexType n = v.length();
257 
258     ASSERT(v.firstIndex()==1);
259     ASSERT(v.length()==n);
260 
261     ASSERT(x.firstIndex()==1);
262     ASSERT(x.length()==n);
263 
264     ASSERT(sgn.firstIndex()==1);
265     ASSERT(sgn.length()==n);
266 
267     ASSERT(iSave.firstIndex()==1);
268     ASSERT(iSave.length()==3);
269 #   endif
270 
271 #   ifdef CHECK_CXXLAPACK
272 //
273 //  Make copies of output arguments
274 //
275     const typename DenseVector<VV>::NoView      v_org     = v; 
276     const typename DenseVector<VX>::NoView      x_org     = x;
277     const typename DenseVector<VSGN>::NoView    sgn_org   = sgn;
278     const EST                                   est_org   = est;
279     const IndexType                             kase_org  = kase; 
280     const typename DenseVector<VSAVE>::NoView   iSave_org = iSave;
281 #   endif
282 
283 //
284 //  Call implementation
285 //
286     lacn2_generic(v, x, sgn, est, kase, iSave);
287 
288 #   ifdef CHECK_CXXLAPACK
289 //
290 //  Make copies of results computed by generic implementation
291 //
292     const typename DenseVector<VV>::NoView      v_generic     = v; 
293     const typename DenseVector<VX>::NoView      x_generic     = x;
294     const typename DenseVector<VSGN>::NoView    sgn_generic   = sgn;
295     const EST                                   est_generic   = est;
296     const IndexType                             kase_generic  = kase; 
297     const typename DenseVector<VSAVE>::NoView   iSave_generic = iSave;
298 
299 //
300 //  restore output arguments
301 //
302     v     = v_org; 
303     x     = x_org;
304     sgn   = sgn_org;
305     est   = est_org;
306     kase  = kase_org; 
307     iSave = iSave_org;
308 
309 //
310 //  Compare generic results with results from the native implementation
311 //
312     lacn2_native(v, x, sgn, est, kase, iSave);
313 
314     bool failed = false;
315     if (! isIdentical(v_generic, v, "v_generic""v")) {
316         std::cerr << "CXXLAPACK: v_generic = " << v_generic << std::endl;
317         std::cerr << "F77LAPACK: v = " << v << std::endl;
318         failed = true;
319     }
320     if (! isIdentical(x_generic, x, "x_generic""x")) {
321         std::cerr << "CXXLAPACK: x_generic = " << x_generic << std::endl;
322         std::cerr << "F77LAPACK: x = " << x << std::endl;
323         failed = true;
324     }
325     if (! isIdentical(sgn_generic, sgn, "sgn_generic""sgn")) {
326         std::cerr << "CXXLAPACK: sgn_generic = " << sgn_generic << std::endl;
327         std::cerr << "F77LAPACK: sgn = " << sgn << std::endl;
328         failed = true;
329     }
330     if (! isIdentical(est_generic, est, "est_generic""est")) {
331         std::cerr << "CXXLAPACK: est_generic = " << est_generic << std::endl;
332         std::cerr << "F77LAPACK: est = " << est << std::endl;
333         failed = true;
334     }
335     if (! isIdentical(kase_generic, kase, "kase_generic""kase")) {
336         std::cerr << "CXXLAPACK: kase_generic = " << kase_generic << std::endl;
337         std::cerr << "F77LAPACK: kase = " << kase << std::endl;
338         failed = true;
339     }
340     if (! isIdentical(iSave_generic, iSave, "iSave_generic""iSave")) {
341         std::cerr << "CXXLAPACK: iSave_generic = "
342                   << iSave_generic << std::endl;
343         std::cerr << "F77LAPACK: iSave = " << iSave << std::endl;
344         failed = true;
345     }
346 
347     if (failed) {
348         std::cerr << "error in: lacn2.tcc" << std::endl;
349         ASSERT(0);
350     } else {
351 //        std::cerr << "passed: lacn2.tcc" << std::endl;
352     }
353 #   endif
354 }
355 
356 //-- forwarding ----------------------------------------------------------------
357 template <typename  VV, typename VX, typename VSGN, typename EST,
358           typename IndexType, typename VSAVE>
359 void
360 lacn2(VV &&v, VX &&x, VSGN &&sgn, EST &&est, IndexType &&kase, VSAVE &&iSave)
361 {
362     lacn2(v, x, sgn, est, kase, iSave);
363 }
364 
365 } } // namespace lapack, flens
366 
367 #endif // FLENS_LAPACK_AUX_LACN2_TCC