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 DGEEQU( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
 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_EQU_TCC
 45 #define FLENS_LAPACK_EIG_EQU_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 template <typename MA, typename VR, typename VC,
 54           typename ROWCOND, typename COLCOND,
 55           typename AMAX>
 56 typename GeMatrix<MA>::IndexType
 57 equ_generic(const GeMatrix<MA>  &A,
 58             DenseVector<VR>     &r,
 59             DenseVector<VC>     &c,
 60             ROWCOND             &rowCond,
 61             COLCOND             &colCond,
 62             AMAX                &amax)
 63 {
 64     using std::abs;
 65     using std::max;
 66     using std::min;
 67 
 68     typedef typename GeMatrix<MA>::ElementType  T;
 69     typedef typename GeMatrix<MA>::IndexType    IndexType;
 70 
 71     const T  Zero(0), One(1);
 72 
 73     const Underscore<IndexType>  _;
 74 
 75     const IndexType m = A.numRows();
 76     const IndexType n = A.numCols();
 77 
 78     IndexType info = 0;
 79 //
 80 //  Quick return if possible
 81 //
 82     if (m==0 || n==0) {
 83         rowCond = One;
 84         colCond = One;
 85         amax = Zero;
 86         return info;
 87     }
 88 //
 89 //  Get machine constants.
 90 //
 91     const T smallNum = lamch<T>(SafeMin);
 92     const T bigNum = One / smallNum;
 93 //
 94 //  Compute row scale factors.
 95 //
 96     r = Zero;
 97 //
 98 //  Find the maximum element in each row.
 99 //
100     for (IndexType j=1; j<=n; ++j) {
101         for (IndexType i=1; i<=m; ++i) {
102             r(i) = max(r(i),abs(A(i,j)));
103         }
104     }
105 //
106 //  Find the maximum and minimum scale factors.
107 //
108     T rcMin = bigNum;
109     T rcMax = Zero;
110     for (IndexType i=1; i<=m; ++i) {
111         rcMax = max(rcMax, r(i));
112         rcMin = min(rcMin, r(i));
113     }
114     amax = rcMax;
115 
116     if (rcMin==Zero) {
117 //
118 //      Find the first zero scale factor and return an error code.
119 //
120         for (IndexType i=1; i<=m; ++i) {
121             if (r(i)==Zero) {
122                 info = i;
123                 return info;
124             }
125         }
126     } else {
127 //
128 //      Invert the scale factors.
129 //
130         for (IndexType i=1; i<=m; ++i) {
131             r(i) = One / min(max(r(i),smallNum), bigNum);
132         }
133 //
134 //      Compute ROWCND = min(R(I)) / max(R(I))
135 //
136         rowCond = max(rcMin,smallNum) / min(rcMax,bigNum);
137     }
138 
139 //  Compute column scale factors
140 //
141     c = Zero;
142 //
143 //  Find the maximum element in each column,
144 //  assuming the row scaling computed above.
145 //
146     for (IndexType j=1; j<=n; ++j) {
147         for (IndexType i=1; i<=m; ++i) {
148             c(j) = max(c(j),abs(A(i,j))*r(i));
149         }
150     }
151 //
152 //  Find the maximum and minimum scale factors.
153 //
154     rcMin = bigNum;
155     rcMax = Zero;
156     for (IndexType j=1; j<=n; ++j) {
157         rcMin = min(rcMin, c(j));
158         rcMax = max(rcMax, c(j));
159     }
160 
161     if (rcMin==Zero) {
162 //
163 //      Find the first zero scale factor and return an error code.
164 //
165         for (IndexType j=1; j<=n; ++j) {
166             if (c(j)==Zero) {
167                 info = m+j;
168                 return info;
169             }
170         }
171     } else {
172 //
173 //      Invert the scale factors.
174 //
175         for (IndexType j=1; j<=n; ++j) {
176             c(j) = One / min(max(c(j),smallNum), bigNum);
177         }
178 //
179 //      Compute COLCND = min(C(J)) / max(C(J))
180 //
181         colCond = max(rcMin,smallNum) / min(rcMax,bigNum);
182     }
183     return info;
184 }
185 
186 
187 //== interface for native lapack ===============================================
188 
189 #ifdef CHECK_CXXLAPACK
190 
191 template <typename MA, typename VR, typename VC,
192           typename ROWCOND, typename COLCOND,
193           typename AMAX>
194 typename GeMatrix<MA>::ElementType
195 equ_native(const GeMatrix<MA>  &A,
196            DenseVector<VR>     &r,
197            DenseVector<VC>     &c,
198            ROWCOND             &rowCond,
199            COLCOND             &colCond,
200            AMAX                &amax)
201 {
202     typedef typename GeMatrix<MA>::ElementType  T;
203 
204     const INTEGER  M = A.numRows();
205     const INTEGER  N = A.numCols();
206     const INTEGER  LDA = A.leadingDimension();
207     INTEGER        INFO;
208 
209     if (IsSame<T,double>::value) {
210         LAPACK_IMPL(dgeequ)(&M,
211                             &N,
212                             A.data(),
213                             &LDA,
214                             r.data(),
215                             c.data(),
216                             &rowCond,
217                             &colCond,
218                             &amax,
219                             &INFO);
220     } else {
221         ASSERT(0);
222     }
223 
224     if (INFO<0) {
225         std::cerr << "dgeequ: INFO = " << INFO << std::endl;
226     }
227     ASSERT(INFO>=0);
228     return INFO;
229 }
230 
231 #endif // CHECK_CXXLAPACK
232 
233 //== public interface ==========================================================
234 template <typename MA, typename VR, typename VC,
235           typename ROWCOND, typename COLCOND,
236           typename AMAX>
237 typename GeMatrix<MA>::IndexType
238 equ(const GeMatrix<MA>  &A,
239     DenseVector<VR>     &r,
240     DenseVector<VC>     &c,
241     ROWCOND             &rowCond,
242     COLCOND             &colCond,
243     AMAX                &amax)
244 {
245     typedef typename GeMatrix<MA>::IndexType  IndexType;
246 //
247 //  Test the input parameters
248 //
249 #   ifndef NDEBUG
250     ASSERT(A.firstRow()==1);
251     ASSERT(A.firstCol()==1);
252 
253     const IndexType m = A.numRows();
254     const IndexType n = A.numCols();
255 
256     ASSERT(r.firstIndex()==1);
257     ASSERT(r.length()==m);
258 
259     ASSERT(c.firstIndex()==1);
260     ASSERT(c.length()==n);
261 #   endif
262 
263 #   ifdef CHECK_CXXLAPACK
264 //
265 //  Make copies of output arguments
266 //
267     typename DenseVector<VR>::NoView  r_org       = r;
268     typename DenseVector<VC>::NoView  c_org       = c;
269     ROWCOND                           rowCond_org = rowCond;
270     COLCOND                           colCond_org = colCond;
271     AMAX                              amax_org    = amax;
272 #   endif
273 
274 //
275 //  Call implementation
276 //
277     const IndexType info = equ_generic(A, r, c, rowCond, colCond, amax);
278 
279 #   ifdef CHECK_CXXLAPACK
280 //
281 //  Make copies of results computed by the generic implementation
282 //
283     typename DenseVector<VR>::NoView  r_generic       = r;
284     typename DenseVector<VC>::NoView  c_generic       = c;
285     ROWCOND                           rowCond_generic = rowCond;
286     COLCOND                           colCond_generic = colCond;
287     AMAX                              amax_generic    = amax;
288 //
289 //  restore output arguments
290 //
291     r = r_org;
292     c = c_org;
293     rowCond = rowCond_org;
294     colCond = colCond_org;
295     amax = amax_org;
296 
297 //
298 //  Compare generic results with results from the native implementation
299 //
300     const IndexType info_ = equ_native(A, r, c, rowCond, colCond, amax);
301 
302     bool failed = false;
303     if (! isIdentical(r_generic, r, "r_generic""r")) {
304         std::cerr << "CXXLAPACK: r_generic = " << r_generic << std::endl;
305         std::cerr << "F77LAPACK: r = " << r << std::endl;
306         failed = true;
307     }
308 
309     if (! isIdentical(c_generic, c, "c_generic""c")) {
310         std::cerr << "CXXLAPACK: c_generic = " << c_generic << std::endl;
311         std::cerr << "F77LAPACK: c = " << c << std::endl;
312         failed = true;
313     }
314 
315     if (! isIdentical(rowCond_generic, rowCond, "rowCond_generic""rowCond")) {
316         std::cerr << "CXXLAPACK: rowCond_generic = "
317                   << rowCond_generic << std::endl;
318         std::cerr << "F77LAPACK: rowCond = "
319                   << rowCond << std::endl;
320         failed = true;
321     }
322 
323     if (! isIdentical(colCond_generic, colCond, "colCond_generic""colCond")) {
324         std::cerr << "CXXLAPACK: colCond_generic = "
325                   << colCond_generic << std::endl;
326         std::cerr << "F77LAPACK: colCond = "
327                   << colCond << std::endl;
328         failed = true;
329     }
330 
331     if (! isIdentical(amax_generic, amax, "amax_generic""amax")) {
332         std::cerr << "CXXLAPACK: amax_generic = "
333                   << amax_generic << std::endl;
334         std::cerr << "F77LAPACK: amax = " << amax << std::endl;
335         failed = true;
336     }
337 
338     if (info!=info_) {
339         std::cerr << "CXXLAPACK: info = " << info << std::endl;
340         std::cerr << "F77LAPACK: info_ = " << info_ << std::endl;
341         failed = true;
342     }
343 
344     if (failed) {
345         std::cerr << "A = " << A << std::endl;
346         std::cerr << "r = " << r << std::endl;
347         std::cerr << "c = " << c << std::endl;
348         ASSERT(0);
349     } else {
350 //        std::cerr << "passed: equ.tcc" << std::endl;
351     }
352 #   endif
353 
354     return info;
355 }
356 
357 
358 //-- forwarding ----------------------------------------------------------------
359 template <typename MA, typename VR, typename VC,
360           typename ROWCOND, typename COLCOND,
361           typename AMAX>
362 typename MA::IndexType
363 equ(const MA  &A,
364     VR        &&r,
365     VC        &&c,
366     ROWCOND   &&rowCond,
367     COLCOND   &&colCond,
368     AMAX      &&amax)
369 {
370     CHECKPOINT_ENTER;
371     const auto info = equ(A, r, c, rowCond, colCond, amax);
372     CHECKPOINT_LEAVE;
373 
374     return info;
375 }
376 
377 } } // namespace lapack, flens
378 
379 #endif // FLENS_LAPACK_EIG_EQU_TCC