1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 2013 Google Inc. All rights reserved.
3 // http://code.google.com/p/ceres-solver/
4 //
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are met:
7 //
8 // * Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright notice,
11 // this list of conditions and the following disclaimer in the documentation
12 // and/or other materials provided with the distribution.
13 // * Neither the name of Google Inc. nor the names of its contributors may be
14 // used to endorse or promote products derived from this software without
15 // specific prior written permission.
16 //
17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 // POSSIBILITY OF SUCH DAMAGE.
28 //
29 // Author: sameeragarwal@google.com (Sameer Agarwal)
30
31 #include "ceres/lapack.h"
32 #include "glog/logging.h"
33
34 // C interface to the LAPACK Cholesky factorization and triangular solve.
35 extern "C" void dpotrf_(char* uplo,
36 int* n,
37 double* a,
38 int* lda,
39 int* info);
40
41 extern "C" void dpotrs_(char* uplo,
42 int* n,
43 int* nrhs,
44 double* a,
45 int* lda,
46 double* b,
47 int* ldb,
48 int* info);
49
50 extern "C" void dgels_(char* uplo,
51 int* m,
52 int* n,
53 int* nrhs,
54 double* a,
55 int* lda,
56 double* b,
57 int* ldb,
58 double* work,
59 int* lwork,
60 int* info);
61
62
63 namespace ceres {
64 namespace internal {
65
SolveInPlaceUsingCholesky(int num_rows,const double * in_lhs,double * rhs_and_solution)66 int LAPACK::SolveInPlaceUsingCholesky(int num_rows,
67 const double* in_lhs,
68 double* rhs_and_solution) {
69 #ifdef CERES_NO_LAPACK
70 LOG(FATAL) << "Ceres was built without a BLAS library.";
71 return -1;
72 #else
73 char uplo = 'L';
74 int n = num_rows;
75 int info = 0;
76 int nrhs = 1;
77 double* lhs = const_cast<double*>(in_lhs);
78
79 dpotrf_(&uplo, &n, lhs, &n, &info);
80 if (info != 0) {
81 LOG(INFO) << "Cholesky factorization (dpotrf) failed: " << info;
82 return info;
83 }
84
85 dpotrs_(&uplo, &n, &nrhs, lhs, &n, rhs_and_solution, &n, &info);
86 if (info != 0) {
87 LOG(INFO) << "Triangular solve (dpotrs) failed: " << info;
88 }
89
90 return info;
91 #endif
92 };
93
EstimateWorkSizeForQR(int num_rows,int num_cols)94 int LAPACK::EstimateWorkSizeForQR(int num_rows, int num_cols) {
95 #ifdef CERES_NO_LAPACK
96 LOG(FATAL) << "Ceres was built without a LAPACK library.";
97 return -1;
98 #else
99 char trans = 'N';
100 int nrhs = 1;
101 int lwork = -1;
102 double work;
103 int info = 0;
104 dgels_(&trans,
105 &num_rows,
106 &num_cols,
107 &nrhs,
108 NULL,
109 &num_rows,
110 NULL,
111 &num_rows,
112 &work,
113 &lwork,
114 &info);
115
116 CHECK_EQ(info, 0);
117 return work;
118 #endif
119 }
120
SolveUsingQR(int num_rows,int num_cols,const double * in_lhs,int work_size,double * work,double * rhs_and_solution)121 int LAPACK::SolveUsingQR(int num_rows,
122 int num_cols,
123 const double* in_lhs,
124 int work_size,
125 double* work,
126 double* rhs_and_solution) {
127 #ifdef CERES_NO_LAPACK
128 LOG(FATAL) << "Ceres was built without a LAPACK library.";
129 return -1;
130 #else
131 char trans = 'N';
132 int m = num_rows;
133 int n = num_cols;
134 int nrhs = 1;
135 int lda = num_rows;
136 int ldb = num_rows;
137 int info = 0;
138 double* lhs = const_cast<double*>(in_lhs);
139
140 dgels_(&trans,
141 &m,
142 &n,
143 &nrhs,
144 lhs,
145 &lda,
146 rhs_and_solution,
147 &ldb,
148 work,
149 &work_size,
150 &info);
151
152 return info;
153 #endif
154 }
155
156 } // namespace internal
157 } // namespace ceres
158