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  */
 36 
 37 #ifndef FLENS_LAPACK_GESV_POTRI_TCC
 38 #define FLENS_LAPACK_GESV_POTRI_TCC 1
 39 
 40 #include <algorithm>
 41 #include <flens/blas/blas.h>
 42 #include <flens/lapack/lapack.h>
 43 
 44 #include <flens/lapack/interface/include/f77lapack.h>
 45 
 46 namespace flens { namespace lapack {
 47 
 48 //== generic lapack implementation =============================================
 49 
 50 template <typename MA>
 51 typename SyMatrix<MA>::IndexType
 52 potri_generic(SyMatrix<MA> &A)
 53 {
 54     typedef typename SyMatrix<MA>::IndexType    IndexType;
 55     const IndexType n = A.dim();
 56     IndexType info = 0;
 57 //
 58 //  Quick return if possible
 59 //
 60     if (n==0) {
 61         return 0;
 62     }
 63 //
 64 //  Invert the triangular Cholesky factor U or L.
 65 //
 66     auto T = A.triangular();
 67 
 68     info = tri(T);
 69     if (info==0) {
 70 //
 71 //      Form inv(U) * inv(U)**T or inv(L)**T * inv(L).
 72 //
 73         lauum(T);
 74     }
 75     return info;
 76 }
 77 
 78 //== interface for native lapack ===============================================
 79 
 80 #ifdef CHECK_CXXLAPACK
 81 
 82 template <typename MA>
 83 typename SyMatrix<MA>::IndexType
 84 potri_native(SyMatrix<MA> &_A)
 85 {
 86     typedef typename SyMatrix<MA>::ElementType  T;
 87 
 88     const char       UPLO = char(_A.upLo());
 89     const INTEGER    N = _A.dim();
 90     T               *A = _A.data();
 91     const INTEGER    LDA = _A.leadingDimension();
 92     INTEGER          INFO;
 93 
 94     if (IsSame<T, DOUBLE>::value) {
 95         LAPACK_IMPL(dpotri)(&UPLO, &N, A, &LDA, &INFO);
 96     } else {
 97         ASSERT(0);
 98     }
 99 
100     return INFO;
101 }
102 
103 #endif // CHECK_CXXLAPACK
104 
105 //== public interface ==========================================================
106 
107 template <typename MA>
108 typename SyMatrix<MA>::IndexType
109 potri(SyMatrix<MA> &A)
110 {
111     typedef typename SyMatrix<MA>::IndexType    IndexType;
112 
113 //
114 //  Test the input parameters
115 //
116 #   ifndef NDEBUG
117     ASSERT(A.firstRow()==1);
118     ASSERT(A.firstCol()==1);
119 #   endif
120 
121 #   ifdef CHECK_CXXLAPACK
122 //
123 //  Make copies of output arguments
124 //
125     typename SyMatrix<MA>::NoView  A_org = A;
126 #   endif
127 
128 //
129 //  Call implementation
130 //
131     const IndexType info = potri_generic(A);
132 
133 #   ifdef CHECK_CXXLAPACK
134 //
135 //  Make copies of generic results
136 //
137     typename SyMatrix<MA>::NoView  A_generic = A;
138 //
139 //  restore output parameters
140 //
141     A = A_org;
142 //
143 //  Compare results
144 //
145     const IndexType _info = potri_native(A);
146 
147     bool failed = false;
148     if (! isIdentical(A_generic, A, "A_generic", "A")) {
149         std::cerr << "A_org =" << A_org << std::endl;
150         std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
151         std::cerr << "F77LAPACK: A = " << A << std::endl;
152         failed = true;
153     }
154 
155     if (! isIdentical(info, _info, " info", "_info")) {
156         std::cerr << "CXXLAPACK:  info = " << info << std::endl;
157         std::cerr << "F77LAPACK: _info = " << _info << std::endl;
158         failed = true;
159     }
160 
161     if (failed) {
162         ASSERT(0);
163     }
164 
165 #   endif
166 
167     return info;
168 }
169 
170 //-- forwarding ----------------------------------------------------------------
171 template <typename MA>
172 typename MA::IndexType
173 potri(MA &&A)
174 {
175     typedef typename MA::IndexType  IndexType;
176 
177     CHECKPOINT_ENTER;
178     const IndexType info =  potri(A);
179     CHECKPOINT_LEAVE;
180 
181     return info;
182 }
183 
184 } } // namespace lapack, flens
185 
186 #endif // FLENS_LAPACK_GESV_POTRI_TCC
Home