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 DPOTF2( UPLO, N, A, LDA, INFO )
 36  *
 37  *  -- LAPACK routine (version 3.3.1) --
 38  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 39  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 40  *  -- April 2011                                                      --
 41  */
 42 
 43 #ifndef FLENS_LAPACK_GESV_POTF2_TCC
 44 #define FLENS_LAPACK_GESV_POTF2_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 
 56 template <typename MA>
 57 typename SyMatrix<MA>::IndexType
 58 potf2_generic(SyMatrix<MA> &A)
 59 {
 60     using std::isnan;
 61     using std::sqrt;
 62 
 63     typedef typename SyMatrix<MA>::ElementType  T;
 64     typedef typename SyMatrix<MA>::IndexType    IndexType;
 65 
 66     const Underscore<IndexType> _;
 67 
 68     const IndexType n = A.dim();
 69     const bool upper = (A.upLo()==Upper);
 70 
 71     const T Zero(0), One(1);
 72 
 73     IndexType info = 0;
 74 //
 75 //  Quick return if possible
 76 //
 77     if (n==0) {
 78         return info;
 79     }
 80     if (upper) {
 81 //
 82 //      Compute the Cholesky factorization A = U**T *U.
 83 //
 84         for (IndexType j=1; j<=n; ++j) {
 85 //
 86 //          Partition matrix
 87 //
 88             const auto range1 = _(1,j-1);
 89             const auto range3 = _(j+1,n);
 90 
 91             const auto a12 = A(range1,j);
 92             const auto A13 = A(range1,range3);
 93             auto a23 = A(j,range3);
 94 //
 95 //          Compute U(J,J) and test for non-positive-definiteness.
 96 //
 97             T a22 = A(j,j) - a12*a12;
 98             if (a22<=Zero || isnan(a22)) {
 99                 A(j,j) = a22;
100                 info = j;
101                 break;
102             }
103             a22 = sqrt(a22);
104             A(j,j) = a22;
105 //
106 //          Compute elements J+1:N of row J.
107 //
108             if (j<n) {
109                 blas::mv(Trans, -One, A13, a12, One, a23);
110                 a23 *= One / a22;
111             }
112         }
113     } else {
114 //
115 //      Compute the Cholesky factorization A = L*L**T.
116 //
117         for (IndexType j=1; j<=n; ++j) {
118 //
119 //          Partition matrix
120 //
121             const auto range1 = _(1,j-1);
122             const auto range3 = _(j+1,n);
123 
124             const auto a21 = A(j,range1);
125             const auto A31 = A(range3,range1);
126             auto a32 = A(range3,j);
127 //
128 //          Compute L(J,J) and test for non-positive-definiteness.
129 //
130             T a22 = A(j,j) - a21*a21;
131             if (a22<=Zero || isnan(a22)) {
132                 A(j,j) = a22;
133                 info = j;
134                 break;
135             }
136             a22 = sqrt(a22);
137             A(j,j) = a22;
138 //
139 //          Compute elements J+1:N of column J.
140 //
141             if (j<n) {
142                 blas::mv(NoTrans, -One, A31, a21, One, a32);
143                 a32 *= One / a22;
144             }
145         }
146     }
147     return info;
148 }
149 
150 //== interface for native lapack ===============================================
151 
152 #ifdef CHECK_CXXLAPACK
153 
154 template <typename MA>
155 typename SyMatrix<MA>::IndexType
156 potf2_native(SyMatrix<MA> &_A)
157 {
158     typedef typename SyMatrix<MA>::ElementType  T;
159 
160     const char       UPLO = char(_A.upLo());
161     const INTEGER    N = _A.dim();
162     T               *A = _A.data();
163     const INTEGER    LDA = _A.leadingDimension();
164     INTEGER          INFO;
165 
166     if (IsSame<T, DOUBLE>::value) {
167         LAPACK_IMPL(dpotf2)(&UPLO, &N, A, &LDA, &INFO);
168     } else {
169         ASSERT(0);
170     }
171 
172     return INFO;
173 }
174 
175 #endif // CHECK_CXXLAPACK
176 
177 //== public interface ==========================================================
178 
179 template <typename MA>
180 typename SyMatrix<MA>::IndexType
181 potf2(SyMatrix<MA> &A)
182 {
183     typedef typename SyMatrix<MA>::IndexType    IndexType;
184 
185 //
186 //  Test the input parameters
187 //
188     ASSERT(A.firstRow()==1);
189     ASSERT(A.firstCol()==1);
190 
191 #   ifdef CHECK_CXXLAPACK
192 //
193 //  Make copies of output arguments
194 //
195     typename SyMatrix<MA>::NoView       A_org      = A;
196 #   endif
197 
198 //
199 //  Call implementation
200 //
201     IndexType info = potf2_generic(A);
202 
203 #   ifdef CHECK_CXXLAPACK
204 //
205 //  Make copies of generic results
206 //
207     typename SyMatrix<MA>::NoView       A_generic      = A;
208 //
209 //  Restore output arguments
210 //
211     A = A_org;
212 
213 //
214 //  Compare results
215 //
216     IndexType _info = potf2_native(A);
217 
218     bool failed = false;
219     if (! isIdentical(A_generic, A, "A_generic", "_A")) {
220         std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
221         std::cerr << "F77LAPACK: A = " << A << std::endl;
222         failed = true;
223     }
224 
225     if (! isIdentical(info, _info, " info", "_info")) {
226         std::cerr << "CXXLAPACK:  info = " << info << std::endl;
227         std::cerr << "F77LAPACK: _info = " << _info << std::endl;
228         failed = true;
229     }
230 
231     if (failed) {
232         ASSERT(0);
233     }
234 
235 #   endif
236 
237     return info;
238 }
239 
240 //-- forwarding ----------------------------------------------------------------
241 template <typename MA>
242 typename MA::IndexType
243 potf2(MA &&A)
244 {
245     typedef typename MA::IndexType  IndexType;
246 
247     CHECKPOINT_ENTER;
248     IndexType info =  potf2(A);
249     CHECKPOINT_LEAVE;
250 
251     return info;
252 }
253 
254 } } // namespace lapack, flens
255 
256 #endif // FLENS_LAPACK_GESV_POTF2_TCC
Home