1 // Copyright (C) 2010  Davis E. King (davis@dlib.net)
2 // License: Boost Software License   See LICENSE.txt for the full license.
3 #ifndef DLIB_LAPACk_GETRF_Hh_
4 #define DLIB_LAPACk_GETRF_Hh_
5 
6 #include "fortran_id.h"
7 #include "../matrix.h"
8 
9 namespace dlib
10 {
11     namespace lapack
12     {
13         namespace binding
14         {
15             extern "C"
16             {
17                 void DLIB_FORTRAN_ID(dgetrf) (const integer *m, const integer *n, double *a,
18                                               const integer *lda, integer *ipiv, integer *info);
19 
20                 void DLIB_FORTRAN_ID(sgetrf) (const integer *m, const integer *n, float *a,
21                                               const integer *lda, integer *ipiv, integer *info);
22 
23             }
24 
getrf(integer m,integer n,double * a,integer lda,integer * ipiv)25             inline int getrf (integer m, integer n, double *a,
26                               integer lda, integer *ipiv)
27             {
28                 integer info = 0;
29                 DLIB_FORTRAN_ID(dgetrf)(&m, &n, a, &lda, ipiv, &info);
30                 return info;
31             }
32 
getrf(integer m,integer n,float * a,integer lda,integer * ipiv)33             inline int getrf (integer m, integer n, float *a,
34                               integer lda, integer *ipiv)
35             {
36                 integer info = 0;
37                 DLIB_FORTRAN_ID(sgetrf)(&m, &n, a, &lda, ipiv, &info);
38                 return info;
39             }
40 
41 
42         }
43 
44     // ------------------------------------------------------------------------------------
45 
46 
47 /*  -- LAPACK routine (version 3.1) -- */
48 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
49 /*     November 2006 */
50 
51 /*     .. Scalar Arguments .. */
52 /*     .. */
53 /*     .. Array Arguments .. */
54 /*     .. */
55 
56 /*  Purpose */
57 /*  ======= */
58 
59 /*  DGETRF computes an LU factorization of a general M-by-N matrix A */
60 /*  using partial pivoting with row interchanges. */
61 
62 /*  The factorization has the form */
63 /*     A = P * L * U */
64 /*  where P is a permutation matrix, L is lower triangular with unit */
65 /*  diagonal elements (lower trapezoidal if m > n), and U is upper */
66 /*  triangular (upper trapezoidal if m < n). */
67 
68 /*  This is the right-looking Level 3 BLAS version of the algorithm. */
69 
70 /*  Arguments */
71 /*  ========= */
72 
73 /*  M       (input) INTEGER */
74 /*          The number of rows of the matrix A.  M >= 0. */
75 
76 /*  N       (input) INTEGER */
77 /*          The number of columns of the matrix A.  N >= 0. */
78 
79 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
80 /*          On entry, the M-by-N matrix to be factored. */
81 /*          On exit, the factors L and U from the factorization */
82 /*          A = P*L*U; the unit diagonal elements of L are not stored. */
83 
84 /*  LDA     (input) INTEGER */
85 /*          The leading dimension of the array A.  LDA >= max(1,M). */
86 
87 /*  IPIV    (output) INTEGER array, dimension (min(M,N)) */
88 /*          The pivot indices; for 1 <= i <= min(M,N), row i of the */
89 /*          matrix was interchanged with row IPIV(i). */
90 
91 /*  INFO    (output) INTEGER */
92 /*          = 0:  successful exit */
93 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
94 /*          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization */
95 /*                has been completed, but the factor U is exactly */
96 /*                singular, and division by zero will occur if it is used */
97 /*                to solve a system of equations. */
98 
99 
100     // ------------------------------------------------------------------------------------
101 
102         template <
103             typename T,
104             long NR1, long NR2,
105             long NC1, long NC2,
106             typename MM,
107             typename layout
108             >
getrf(matrix<T,NR1,NC1,MM,column_major_layout> & a,matrix<integer,NR2,NC2,MM,layout> & ipiv)109         int getrf (
110             matrix<T,NR1,NC1,MM,column_major_layout>& a,
111             matrix<integer,NR2,NC2,MM,layout>& ipiv
112         )
113         {
114             const long m = a.nr();
115             const long n = a.nc();
116 
117             ipiv.set_size(std::min(m,n), 1);
118 
119             // compute the actual decomposition
120             return binding::getrf(m, n, &a(0,0), a.nr(), &ipiv(0,0));
121         }
122 
123     // ------------------------------------------------------------------------------------
124 
125     }
126 
127 }
128 
129 // ----------------------------------------------------------------------------------------
130 
131 #endif // DLIB_LAPACk_GETRF_Hh_
132 
133