1 /*
  2  *   Copyright (c) 2009, 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_BLAS_LEVEL1_COPY_TCC
 34 #define FLENS_BLAS_LEVEL1_COPY_TCC 1
 35 
 36 #include <flens/blas/blas.h>
 37 #include <flens/typedefs.h>
 38 
 39 #ifdef FLENS_DEBUG_CLOSURES
 40 #   include <flens/blas/blaslogon.h>
 41 #else
 42 #   include <flens/blas/blaslogoff.h>
 43 #endif
 44 
 45 
 46 namespace flens { namespace blas {
 47 
 48 //-- copy
 49 template <typename VX, typename VY>
 50 void
 51 copy(const DenseVector<VX> &x, DenseVector<VY> &y)
 52 {
 53     FLENS_BLASLOG_SETTAG("--> ");
 54     FLENS_BLASLOG_BEGIN_COPY(x, y);
 55 
 56 //
 57 //  Resize left hand size if needed.  This is *usually* only alloweded
 58 //  when the left hand side is an empty vector (such that it is no actual
 59 //  resizing but rather an initialization).
 60 //
 61     if (y.length()!=x.length()) {
 62 #       ifndef FLENS_DEBUG_CLOSURES
 63         ASSERT(y.length()==0);
 64 #       else
 65         if (y.length()!=0) {
 66             FLENS_BLASLOG_RESIZE_VECTOR(y, x.length());
 67         }
 68 #       endif
 69         y.resize(x);
 70     }
 71     y.changeIndexBase(x.firstIndex());
 72 
 73 #   ifdef HAVE_CXXBLAS_COPY
 74     cxxblas::copy(x.length(), x.data(), x.stride(), y.data(), y.stride());
 75 #   else
 76     ASSERT(0);
 77 #   endif
 78 
 79     FLENS_BLASLOG_END;
 80     FLENS_BLASLOG_UNSETTAG;
 81 }
 82 
 83 //-- gecopy
 84 template <typename MA, typename MB>
 85 void
 86 copy(Transpose trans, const GeMatrix<MA> &A, GeMatrix<MB> &B)
 87 {
 88 //
 89 //  check if this is an inplace transpose of A
 90 //
 91     if (trans==Trans || trans==ConjTrans) {
 92         if (DEBUGCLOSURE::identical(A, B)) {
 93 #           ifndef FLENS_DEBUG_CLOSURES
 94 //
 95 //          temporaries are not allowed
 96 //
 97             ASSERT(A.numRows()==A.numCols());
 98             gecotr(MB::order, trans, B.numRows(), B.numCols(),
 99                    B.data(), B.leadingDimension());
100 #           else
101 //
102 //          temporaries are allowed: check if this requires a temporary
103 //
104             if (A.numRows()!=A.numCols()) {
105                 typename Result<GeMatrix<MA> >::Type _A = A;
106                 FLENS_BLASLOG_TMP_ADD(_A);
107 
108                 copy(trans, _A, B);
109 
110                 FLENS_BLASLOG_TMP_REMOVE(_A, A);
111                 return;
112             } else {
113 //
114 //              otherwise perform inplace transpose
115 //
116                 gecotr(MB::order, trans, B.numRows(), B.numCols(),
117                        B.data(), B.leadingDimension());
118                 return;
119             }
120 #           endif
121         }
122     }
123 
124 //
125 //  Resize left hand size if needed.  This is *usually* only alloweded
126 //  when the left hand side is an empty matrix (such that it is no actual
127 //  resizing but rather an initialization).
128 //
129     if (trans==NoTrans) {
130         if ((A.numRows()!=B.numRows()) || (A.numCols()!=B.numCols())) {
131 #           ifndef FLENS_DEBUG_CLOSURES
132             ASSERT(B.numRows()==0 || B.numCols()==0);
133 #           else
134             if (B.numRows()!=0 || B.numCols()!=0) {
135                 FLENS_BLASLOG_RESIZE_MATRIX(B, A.numRows(), A.numCols());
136             }
137 #           endif
138             B.resize(A);
139         }
140     } else {
141         if ((A.numRows()!=B.numCols())  || (A.numCols()!=B.numRows())) {
142 #           ifndef FLENS_DEBUG_CLOSURES
143             ASSERT(B.numRows()==0 && A.numCols()==0);
144 #           else
145             if (B.numRows()!=0 || B.numCols()!=0) {
146                 FLENS_BLASLOG_RESIZE_MATRIX(B, A.numCols(), A.numRows());
147             }
148 #           endif
149             B.resize(A.numCols(), A.numRows(),
150                      A.firstCol(), A.firstRow());
151         }
152     }
153 
154     trans = (MA::order==MB::order)
155           ? Transpose(trans ^ NoTrans)
156           : Transpose(trans ^ Trans);
157 
158     FLENS_BLASLOG_SETTAG("--> ");
159     FLENS_BLASLOG_BEGIN_MCOPY(trans, A, B);
160 
161 #   ifdef HAVE_CXXBLAS_GECOPY
162     cxxblas::gecopy(MB::order, trans,
163                     B.numRows(), B.numCols(),
164                     A.data(), A.leadingDimension(),
165                     B.data(), B.leadingDimension());
166 #   else
167     ASSERT(0);
168 #   endif
169 
170     FLENS_BLASLOG_END;
171     FLENS_BLASLOG_UNSETTAG;
172 }
173 
174 //-- trcopy
175 template <typename MA, typename MB>
176 void
177 copy(Transpose trans, const TrMatrix<MA> &A, TrMatrix<MB> &B)
178 {
179 //
180 //  Resize left hand size if needed.  This is *usually* only alloweded
181 //  when the left hand side is an empty matrix (such that it is no actual
182 //  resizing but rather an initialization).
183 //
184     if (trans==NoTrans) {
185         if ((A.numRows()!=B.numRows()) || (A.numCols()!=B.numCols())) {
186 #           ifndef FLENS_DEBUG_CLOSURES
187             ASSERT(B.numRows()==0 || B.numCols()==0);
188 #           else
189             if (B.numRows()!=0 && B.numCols()!=0) {
190                 FLENS_BLASLOG_RESIZE_MATRIX(B, A.numRows(), A.numCols());
191             }
192 #           endif
193             B.resize(A);
194         }
195     } else {
196         if ((A.numRows()!=B.numCols()) || (A.numCols()!=B.numRows())) {
197 #           ifndef FLENS_DEBUG_CLOSURES
198             ASSERT(B.numRows()==0 || B.numCols()==0);
199 #           else
200             if (B.numRows()!=0 && B.numCols()!=0) {
201                 FLENS_BLASLOG_RESIZE_MATRIX(B, A.numCols(), A.numRows());
202             }
203 #           endif
204             B.resize(A.numCols(), A.numRows(),
205                      A.firstCol(), A.firstRow());
206         }
207     }
208 
209 #   ifndef NDEBUG
210     if (trans==NoTrans || trans==Conj) {
211         ASSERT(A.upLo()==B.upLo());
212     } else {
213         ASSERT(A.upLo()!=B.upLo());
214     }
215 #   endif
216 
217     // TODO: make this assertion unnecessary
218     ASSERT(A.order()==B.order());
219     ASSERT(A.diag()==B.diag());
220 
221     FLENS_BLASLOG_SETTAG("--> ");
222     FLENS_BLASLOG_BEGIN_MCOPY(trans, A, B);
223 
224 #   ifdef HAVE_CXXBLAS_TRCOPY
225     cxxblas::trcopy(B.order(), B.upLo(), trans, B.diag(),
226                     B.numRows(), B.numCols(),
227                     A.data(), A.leadingDimension(),
228                     B.data(), B.leadingDimension());
229 #   else
230     ASSERT(0);
231 #   endif
232 
233     FLENS_BLASLOG_END;
234     FLENS_BLASLOG_UNSETTAG;
235 }
236 
237 //-- sycopy
238 template <typename MA, typename MB>
239 void
240 copy(const SyMatrix<MA> &A, SyMatrix<MB> &B)
241 {
242 //
243 //  Resize left hand size if needed.  This is *usually* only alloweded
244 //  when the left hand side is an empty matrix (such that it is no actual
245 //  resizing but rather an initialization).
246 //
247     if (B.dim()!=A.dim()) {
248 #       ifndef FLENS_DEBUG_CLOSURES
249         ASSERT(B.dim()==0);
250 #       else
251         if (B.dim()!=0) {
252             FLENS_BLASLOG_RESIZE_MATRIX(B, A.dim(), A.dim());
253         }
254 #       endif
255         B.resize(A);
256         B.upLo() = A.upLo();
257     }
258 
259     ASSERT(A.upLo()==B.upLo());
260     // TODO: make this assertion unnecessary
261     ASSERT(A.order()==B.order());
262 
263     FLENS_BLASLOG_SETTAG("--> ");
264     FLENS_BLASLOG_BEGIN_COPY(A, B);
265 
266 #   ifdef HAVE_CXXBLAS_TRCOPY
267     cxxblas::trcopy(B.order(), B.upLo(), NoTrans, NonUnit, B.dim(), B.dim(),
268                     A.data(), A.leadingDimension(),
269                     B.data(), B.leadingDimension());
270 #   else
271     ASSERT(0);
272 #   endif
273     FLENS_BLASLOG_END;
274     FLENS_BLASLOG_UNSETTAG;
275 }
276 
277 //-- extensions ----------------------------------------------------------------
278 
279 //-- copy: TrMatrix -> GeMatrix
280 template <typename MA, typename MB>
281 void
282 copy(Transpose trans, const TrMatrix<MA> &A, GeMatrix<MB> &B)
283 {
284     typename GeMatrix<MB>::ElementType  Zero(0);
285 
286     if (trans==NoTrans) {
287         if (A.numRows()!=B.numRows() && A.numCols()!=B.numCols()) {
288 #           ifndef FLENS_DEBUG_CLOSURES
289             ASSERT(B.numRows()==0 || B.numCols()==0);
290 #           else
291             if (B.numRows()!=0 && B.numCols()!=0) {
292                 FLENS_BLASLOG_RESIZE_MATRIX(B, A.numRows(), A.numCols());
293             }
294 #           endif
295             B.resize(A.numRows(), A.numCols(),
296                      A.firstRow(), A.firstCol());
297         }
298     } else {
299         if (A.numRows()!=B.numCols() && A.numCols()!=B.numRows()) {
300 #           ifndef FLENS_DEBUG_CLOSURES
301             ASSERT(B.numRows()==0 || B.numCols()==0);
302 #           else
303             if (B.numRows()!=0 && B.numCols()!=0) {
304                 FLENS_BLASLOG_RESIZE_MATRIX(B, A.numCols(), A.numRows());
305             }
306 #           endif
307             B.resize(A.numCols(), A.numRows(),
308                      A.firstCol(), A.firstRow());
309          }
310     }
311 
312     if (trans==NoTrans) {
313         if (A.upLo()==Upper) {
314             B.upper() = A;
315             B.strictLower() = Zero;
316         } else {
317             B.lower() = A;
318             B.strictUpper() = Zero;
319         }
320     } else if (trans==Trans) {
321         if (A.upLo()==Upper) {
322             B.lower() = transpose(A);
323             B.strictUpper() = Zero;
324         } else {
325             B.upper() = transpose(A);
326             B.strictLower() = Zero;
327         }
328     } else {
329         ASSERT(0);
330     }
331 }
332 
333 //-- copy: SyMatrix -> GeMatrix
334 template <typename MA, typename MB>
335 void
336 copy(const SyMatrix<MA> &A, GeMatrix<MB> &B)
337 {
338     if (A.numRows()!=B.numRows() && A.numCols()!=B.numCols()) {
339 #       ifndef FLENS_DEBUG_CLOSURES
340         ASSERT(B.numRows()==0 || B.numCols()==0);
341 #       else
342         if (B.numRows()!=0 && B.numCols()!=0) {
343             FLENS_BLASLOG_RESIZE_MATRIX(B, A.numRows(), A.numCols());
344         }
345 #       endif
346         B.resize(A.numRows(), A.numCols(),
347                  A.firstRow(), A.firstCol());
348     }
349 
350     if (A.upLo()==Upper) {
351         B.upper() = A.general().upper();
352         B.strictLower() = transpose(A.general().strictUpper());
353     } else {
354         B.lower() = A.general().lower();
355         B.strictUpper() = transpose(A.general().strictLower());
356     }
357 }
358 
359 
360 } } // namespace blas, flens
361 
362 #endif // FLENS_BLAS_LEVEL1_COPY_TCC