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 #ifndef FLENS_LAPACK_DEBUG_ISIDENTICAL_TCC
 34 #define FLENS_LAPACK_DEBUG_ISIDENTICAL_TCC 1
 35 
 36 namespace flens { namespace lapack {
 37 
 38 template <typename X, typename Y>
 39 bool
 40 isDifferent(const X &x, const Y &y)
 41 {
 42     using std::isnan;
 43     
 44     if (isnan(x) && isnan(y)) {
 45         return false;
 46     }
 47     return x!=y;
 48 }
 49 
 50 template <typename X, typename Y>
 51 bool
 52 isIdentical(const X &x, const Y &y, const char *xName, const char *yName)
 53 {
 54     if (isDifferent(x,y)) {
 55         std::cerr.precision(150);
 56         std::cerr << xName << " = " << x
 57                   << std::endl
 58                   << yName << " = " << y
 59                   << std::endl
 60                   << xName << " - " << yName << " = "
 61                   << x - y
 62                   << std::endl;
 63             return false;
 64     }
 65     return true;
 66 }
 67 
 68 template <typename VX, typename VY>
 69 bool
 70 isIdentical(const DenseVector<VX> &x, const DenseVector<VY> &y,
 71             const char *xName, const char *yName)
 72 {
 73     typedef typename DenseVector<VX>::IndexType IndexType;
 74 
 75     if (x.length()!=y.length()) {
 76         std::cerr << xName << ".length() = " << x.length() << ", "
 77                   << yName << ".length() = " << y.length()
 78                   << std::endl;
 79         return false;
 80     }
 81     if (x.firstIndex()!=y.firstIndex()) {
 82         std::cerr << xName << ".firstIndex() = " << x.firstIndex() << ", "
 83                   << yName << ".firstIndex() = " << y.firstIndex()
 84                   << std::endl;
 85         return false;
 86     }
 87     if (x.endIndex()!=y.endIndex()) {
 88         std::cerr << xName << ".endIndex() = " << x.endIndex() << ", "
 89                   << yName << ".endIndex() = " << y.endIndex()
 90                   << std::endl;
 91         return false;
 92     }
 93     if (x.inc()!=y.inc()) {
 94         std::cerr << xName << ".inc() = " << x.inc() << ", "
 95                   << yName << ".inc() = " << y.inc()
 96                   << std::endl;
 97         return false;
 98     }
 99 
100     for (IndexType i=x.firstIndex(); i!=x.endIndex(); i+=x.inc()) {
101         if (isDifferent(x(i), y(i))) {
102             std::cerr.precision(150);
103             std::cerr << xName << "(" << i << ") = " << x(i)
104                       << std::endl
105                       << yName << "(" << i << ") = " << y(i)
106                       << std::endl
107                       << xName << "(" << i << ") - "
108                       << yName << "(" << i << ") = "
109                       << x(i) - y(i)
110                       << std::endl;
111             return false;
112         }
113     }
114     return true;
115 }
116 
117 template <typename MA, typename MB>
118 bool
119 isIdentical(const GeMatrix<MA> &A, const GeMatrix<MB> &B,
120             const char *AName, const char *BName)
121 {
122     typedef typename GeMatrix<MA>::IndexType IndexType;
123 
124     if (A.numRows()*A.numCols()==0 && B.numRows()*B.numCols()==0) {
125         return true;
126     }
127 
128     if (A.numRows()!=B.numRows()) {
129         std::cerr << AName << ".numRows() = " << A.numRows() << ", "
130                   << BName << ".numRows() = " << B.numRows()
131                   << std::endl;
132         return false;
133     }
134     if (A.numCols()!=B.numCols()) {
135         std::cerr << AName << ".numCols() = " << A.numCols() << ", "
136                   << BName << ".numCols() = " << B.numCols()
137                   << std::endl;
138         return false;
139     }
140     if (A.firstRow()!=B.firstRow()) {
141         std::cerr << AName << ".firstRow() = " << A.firstRow() << ", "
142                   << BName << ".firstRow() = " << B.firstRow()
143                   << std::endl;
144         return false;
145     }
146     if (A.firstCol()!=B.firstCol()) {
147         std::cerr << AName << ".firstCol() = " << A.firstCol() << ", "
148                   << BName << ".firstCol() = " << B.firstCol()
149                   << std::endl;
150         return false;
151     }
152     if (A.lastRow()!=B.lastRow()) {
153         std::cerr << AName << ".lastRow() = " << A.lastRow() << ", "
154                   << BName << ".lastRow() = " << B.lastRow()
155                   << std::endl;
156         return false;
157     }
158     if (A.lastCol()!=B.lastCol()) {
159         std::cerr << AName << ".lastCol() = " << A.lastCol() << ", "
160                   << BName << ".lastCol() = " << B.lastCol()
161                   << std::endl;
162         return false;
163     }
164 
165     for (IndexType i=A.firstRow(); i<=A.lastRow(); ++i) {
166         for (IndexType j=A.firstCol(); j<=A.lastCol(); ++j) {
167             if (isDifferent(A(i,j), B(i,j))) {
168                 std::cerr.precision(150);
169                 std::cerr << AName << "(" << i << ", " << j << ") = " << A(i,j)
170                           << std::endl
171                           << BName << "(" << i << ", " << j << ") = " << B(i,j)
172                           << std::endl
173                           << AName << "(" << i << ", " << j << ") - "
174                           << BName << "(" << i << ", " << j << ") = "
175                           << A(i,j) - B(i,j)
176                           << std::endl;
177                 return false;
178             }
179         }
180     }
181     return true;
182 }
183 
184 template <typename MA, typename MB>
185 bool
186 isIdentical(const TrMatrix<MA> &A, const TrMatrix<MB> &B,
187             const char *AName, const char *BName)
188 {
189     typedef typename TrMatrix<MA>::IndexType IndexType;
190 
191     if (A.dim()==0 && B.dim()==0) {
192         return true;
193     }
194 
195     if (A.dim()!=B.dim()) {
196         std::cerr << AName << ".dim() = " << A.dim() << ", "
197                   << BName << ".dim() = " << B.dim()
198                   << std::endl;
199         return false;
200     }
201     if (A.firstRow()!=B.firstRow()) {
202         std::cerr << AName << ".firstRow() = " << A.firstRow() << ", "
203                   << BName << ".firstRow() = " << B.firstRow()
204                   << std::endl;
205         return false;
206     }
207     if (A.firstCol()!=B.firstCol()) {
208         std::cerr << AName << ".firstCol() = " << A.firstCol() << ", "
209                   << BName << ".firstCol() = " << B.firstCol()
210                   << std::endl;
211         return false;
212     }
213     if (A.lastRow()!=B.lastRow()) {
214         std::cerr << AName << ".lastRow() = " << A.lastRow() << ", "
215                   << BName << ".lastRow() = " << B.lastRow()
216                   << std::endl;
217         return false;
218     }
219     if (A.lastCol()!=B.lastCol()) {
220         std::cerr << AName << ".lastCol() = " << A.lastCol() << ", "
221                   << BName << ".lastCol() = " << B.lastCol()
222                   << std::endl;
223         return false;
224     }
225 
226     if (A.upLo()!=B.upLo()) {
227         std::cerr << AName << ".upLo() = " << A.upLo() << ", "
228                   << BName << ".upLo() = " << B.upLo()
229                   << std::endl;
230         return false;
231     }
232 
233     bool failed = false;
234     bool isUpper = (A.upLo()==Upper);
235     bool isUnit = (A.diag()==Unit);
236 
237     for (IndexType i=A.firstRow(); i<=A.lastRow(); ++i) {
238         for (IndexType j=A.firstCol(); j<=A.lastCol(); ++j) {
239             if (isUnit && i==j) {
240                 continue;
241             }
242             if (isUpper && j>=i) {
243                 failed = isDifferent(A(i,j), B(i,j));
244             }
245             if (!isUpper && j<=i) {
246                 failed = isDifferent(A(i,j), B(i,j));
247             }
248             if (failed) {
249                 std::cerr.precision(150);
250                 std::cerr << AName << "(" << i << ", " << j << ") = "
251                           << A(i,j) << std::endl
252                           << BName << "(" << i << ", " << j << ") = "
253                           << B(i,j) << std::endl
254                           << AName << "(" << i << ", " << j << ") - "
255                           << BName << "(" << i << ", " << j << ") = "
256                           << A(i,j) - B(i,j)
257                           << std::endl;
258                 return false;
259             }
260         }
261     }
262     return true;
263 }
264 
265 template <typename MA, typename MB>
266 bool
267 isIdentical(const SyMatrix<MA> &A, const SyMatrix<MB> &B,
268             const char *AName, const char *BName)
269 {
270     typedef typename SyMatrix<MA>::IndexType IndexType;
271 
272     if (A.dim()==0 && B.dim()==0) {
273         return true;
274     }
275 
276     if (A.dim()!=B.dim()) {
277         std::cerr << AName << ".dim() = " << A.dim() << ", "
278                   << BName << ".dim() = " << B.dim()
279                   << std::endl;
280         return false;
281     }
282     if (A.firstRow()!=B.firstRow()) {
283         std::cerr << AName << ".firstRow() = " << A.firstRow() << ", "
284                   << BName << ".firstRow() = " << B.firstRow()
285                   << std::endl;
286         return false;
287     }
288     if (A.firstCol()!=B.firstCol()) {
289         std::cerr << AName << ".firstCol() = " << A.firstCol() << ", "
290                   << BName << ".firstCol() = " << B.firstCol()
291                   << std::endl;
292         return false;
293     }
294     if (A.lastRow()!=B.lastRow()) {
295         std::cerr << AName << ".lastRow() = " << A.lastRow() << ", "
296                   << BName << ".lastRow() = " << B.lastRow()
297                   << std::endl;
298         return false;
299     }
300     if (A.lastCol()!=B.lastCol()) {
301         std::cerr << AName << ".lastCol() = " << A.lastCol() << ", "
302                   << BName << ".lastCol() = " << B.lastCol()
303                   << std::endl;
304         return false;
305     }
306 
307     if (A.upLo()!=B.upLo()) {
308         std::cerr << AName << ".upLo() = " << A.upLo() << ", "
309                   << BName << ".upLo() = " << B.upLo()
310                   << std::endl;
311         return false;
312     }
313 
314     bool failed = false;
315     bool isUpper = (A.upLo()==Upper);
316 
317     for (IndexType i=A.firstRow(); i<=A.lastRow(); ++i) {
318         for (IndexType j=A.firstCol(); j<=A.lastCol(); ++j) {
319             if (isUpper && j>=i) {
320                 failed = isDifferent(A(i,j), B(i,j));
321             }
322             if (!isUpper && j<=i) {
323                 failed = isDifferent(A(i,j), B(i,j));
324             }
325             if (failed) {
326                 std::cerr.precision(150);
327                 std::cerr << AName << "(" << i << ", " << j << ") = "
328                           << A(i,j) << std::endl
329                           << BName << "(" << i << ", " << j << ") = "
330                           << B(i,j) << std::endl
331                           << AName << "(" << i << ", " << j << ") - "
332                           << BName << "(" << i << ", " << j << ") = "
333                           << A(i,j) - B(i,j)
334                           << std::endl;
335                 return false;
336             }
337         }
338     }
339     return true;
340 }
341 
342 
343 } } // namespace lapack, flens
344 
345 #endif // FLENS_LAPACK_DEBUG_ISIDENTICAL_TCC