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 DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
 36  *
 37  *  -- LAPACK 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 #ifndef FLENS_LAPACK_GESV_TI2_TCC
 44 #define FLENS_LAPACK_GESV_TI2_TCC 1
 45 
 46 #include <algorithm>
 47 #include <flens/blas/blas.h>
 48 #include <flens/lapack/lapack.h>
 49 
 50 #include <flens/lapack/interface/include/f77lapack.h>
 51 
 52 namespace flens { namespace lapack {
 53 
 54 //== generic lapack implementation =============================================
 55 template <typename MA>
 56 typename TrMatrix<MA>::IndexType
 57 ti2_generic(TrMatrix<MA> &A)
 58 {
 59     typedef typename TrMatrix<MA>::ElementType ElementType;
 60     typedef typename TrMatrix<MA>::IndexType   IndexType;
 61 
 62     const Underscore<IndexType>  _;
 63 
 64     const bool upper  = (A.upLo()==Upper);
 65     const bool noUnit = (A.diag()==NonUnit);
 66     const IndexType n = A.dim();
 67 
 68     const ElementType  One(1);
 69 
 70     IndexType info = 0;
 71 
 72     if (upper) {
 73 //
 74 //      Compute inverse of upper triangular matrix.
 75 //
 76         for (IndexType j=1; j<=n; ++j) {
 77             ElementType Ajj;
 78             if (noUnit) {
 79                 A(j,j) = One / A(j,j);
 80                 Ajj = -A(j,j);
 81             } else {
 82                 Ajj = -One;
 83             }
 84 //
 85 //          Compute elements 1:j-1 of j-th column.
 86 //
 87             const auto range = _(1,j-1);
 88 
 89             const auto U11 = (noUnit) ? A(range,range).upper()
 90                                       : A(range,range).upperUnit();
 91             auto u12 = A(range,j);
 92 
 93             blas::mv(NoTrans, U11, u12);
 94             u12 *= Ajj;
 95         }
 96     } else {
 97 //
 98 //      Compute inverse of lower triangular matrix.
 99 //
100         for (IndexType j=n; j>=1; --j) {
101             ElementType Ajj;
102             if (noUnit) {
103                 A(j,j) = One / A(j,j);
104                 Ajj = -A(j,j);
105             } else {
106                 Ajj = -One;
107             }
108             if (j<n) {
109 //
110 //              Compute elements j+1:n of j-th column.
111 //
112                 const auto range = _(j+1,n);
113 
114                 const auto L22 = (noUnit) ? A(range,range).lower()
115                                           : A(range,range).lowerUnit();
116                 auto l21 = A(range,j);
117 
118                 blas::mv(NoTrans, L22, l21);
119                 l21 *= Ajj;
120             }
121         }
122     }
123 
124     return info;
125 }
126 
127 //== interface for native lapack ===============================================
128 
129 #ifdef CHECK_CXXLAPACK
130 
131 template <typename MA>
132 typename TrMatrix<MA>::IndexType
133 ti2_native(TrMatrix<MA> &A)
134 {
135     typedef typename GeMatrix<MA>::ElementType  ElementType;
136 
137     const char       UPLO   = getF77BlasChar(A.upLo());
138     const char       DIAG   = getF77BlasChar(A.diag());
139     const INTEGER    N      = A.dim();
140     const INTEGER    LDA    = A.leadingDimension();
141     INTEGER          INFO;
142 
143     if (IsSame<ElementType, DOUBLE>::value) {
144         LAPACK_IMPL(dtrti2)(&UPLO,
145                             &DIAG,
146                             &N,
147                             A.data(),
148                             &LDA,
149                             &INFO);
150     } else {
151         ASSERT(0);
152     }
153     ASSERT(INFO>=0);
154     return INFO;
155 }
156 
157 #endif // CHECK_CXXLAPACK
158 
159 //== public interface ==========================================================
160 
161 template <typename MA>
162 typename TrMatrix<MA>::IndexType
163 ti2(TrMatrix<MA> &A)
164 {
165     typedef typename GeMatrix<MA>::IndexType IndexType;
166 
167 //
168 //  Test the input parameters
169 //
170 #   ifndef NDEBUG
171     ASSERT(A.firstRow()==1);
172     ASSERT(A.firstCol()==1);
173 #   endif
174 //
175 //  Make copies of output arguments
176 //
177 #   ifdef CHECK_CXXLAPACK
178     typename TrMatrix<MA>::NoView   A_org = A;
179 #   endif
180 
181 //
182 //  Call implementation
183 //
184     const IndexType info = ti2_generic(A);
185 
186 //
187 //  Compare results
188 //
189 #   ifdef CHECK_CXXLAPACK
190     typename TrMatrix<MA>::NoView   A_generic = A;
191 
192     A = A_org;
193 
194     const IndexType _info = ti2_native(A);
195 
196     bool failed = false;
197     if (! isIdentical(A_generic, A, "A_generic", "A")) {
198         std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
199         std::cerr << "F77LAPACK: A = " << A << std::endl;
200         failed = true;
201     }
202 
203     if (! isIdentical(info, _info, " info", "_info")) {
204         std::cerr << "CXXLAPACK:  info = " << info << std::endl;
205         std::cerr << "F77LAPACK: _info = " << _info << std::endl;
206         failed = true;
207     }
208 
209     if (failed) {
210         ASSERT(0);
211     }
212 #   endif
213 
214     return info;
215 }
216 
217 //-- forwarding ----------------------------------------------------------------
218 
219 template <typename MA>
220 typename MA::IndexType
221 ti2(MA &&A)
222 {
223     typedef typename MA::IndexType  IndexType;
224 
225     CHECKPOINT_ENTER;
226     IndexType info = ti2(A);
227     CHECKPOINT_LEAVE;
228 
229     return info;
230 }
231 
232 } } // namespace lapack, flens
233 
234 #endif // FLENS_LAPACK_GESV_TI2_TCC
Home