1 /*
2  *    This file is part of CasADi.
3  *
4  *    CasADi -- A symbolic framework for dynamic optimization.
5  *    Copyright (C) 2010-2014 Joel Andersson, Joris Gillis, Moritz Diehl,
6  *                            K.U. Leuven. All rights reserved.
7  *    Copyright (C) 2011-2014 Greg Horn
8  *
9  *    CasADi is free software; you can redistribute it and/or
10  *    modify it under the terms of the GNU Lesser General Public
11  *    License as published by the Free Software Foundation; either
12  *    version 3 of the License, or (at your option) any later version.
13  *
14  *    CasADi is distributed in the hope that it will be useful,
15  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  *    Lesser General Public License for more details.
18  *
19  *    You should have received a copy of the GNU Lesser General Public
20  *    License along with CasADi; if not, write to the Free Software
21  *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22  *
23  */
24 
25 
26 #ifndef CASADI_LAPACK_QR_HPP
27 #define CASADI_LAPACK_QR_HPP
28 
29 #include "casadi/core/linsol_internal.hpp"
30 #include <casadi/interfaces/lapack/casadi_linsol_lapackqr_export.h>
31 
32 extern "C" {
33   /// QR-factorize dense matrix (lapack)
34   void dgeqrf_(int *m, int *n, double *a, int *lda, double *tau,
35                double *work, int *lwork, int *info);
36 
37   /// Multiply right hand side with Q-transpose (lapack)
38   void dormqr_(char *side, char *trans, int *n, int *m, int *k, double *a,
39                int *lda, double *tau, double *c, int *ldc,
40                double *work, int *lwork, int *info);
41 
42   /// Solve upper triangular system (lapack)
43   void dtrsm_(char *side, char *uplo, char *transa, char *diag, int *m, int *n,
44                          double *alpha, double *a, int *lda, double *b, int *ldb);
45 }
46 
47 /** \defgroup plugin_Linsol_lapackqr
48 *
49 * This class solves the linear system <tt>A.x=b</tt> by making an QR factorization of A: \n
50 * <tt>A = Q.R</tt>, with Q orthogonal and R upper triangular
51 */
52 
53 /** \pluginsection{Linsol,lapackqr} */
54 
55 /// \cond INTERNAL
56 namespace casadi {
57   struct CASADI_LINSOL_LAPACKQR_EXPORT LapackQrMemory : public LinsolMemory {
58     // Matrix
59     std::vector<double> mat;
60 
61     // The scalar factors of the elementary reflectors
62     std::vector<double> tau;
63 
64     // qr work array
65     std::vector<double> work;
66   };
67 
68   /** \brief  \pluginbrief{Linsol,lapackqr}
69    *
70    @copydoc Linsol_doc
71    @copydoc plugin_Linsol_lapackqr
72    *
73    */
74   class CASADI_LINSOL_LAPACKQR_EXPORT LapackQr : public LinsolInternal {
75   public:
76     // Create a linear solver given a sparsity pattern and a number of right hand sides
77     LapackQr(const std::string& name, const Sparsity& sp);
78 
79     /** \brief  Create a new Linsol */
creator(const std::string & name,const Sparsity & sp)80     static LinsolInternal* creator(const std::string& name, const Sparsity& sp) {
81       return new LapackQr(name, sp);
82     }
83 
84     // Destructor
85     ~LapackQr() override;
86 
87     // Initialize the solver
88     void init(const Dict& opts) override;
89 
90     ///@{
91     /** \brief Options */
92     static const Options options_;
get_options() const93     const Options& get_options() const override { return options_;}
94     ///@}
95 
96     /** \brief Create memory block */
alloc_mem() const97     void* alloc_mem() const override { return new LapackQrMemory();}
98 
99     /** \brief Initalize memory block */
100     int init_mem(void* mem) const override;
101 
102     /** \brief Free memory block */
free_mem(void * mem) const103     void free_mem(void *mem) const override { delete static_cast<LapackQrMemory*>(mem);}
104 
105     // Factorize the linear system
106     int nfact(void* mem, const double* A) const override;
107 
108     // Solve the linear system
109     int solve_batch(void* mem, const double* A, double* x, casadi_int nrhs, bool tr) const;
110 
111     // Solve the linear system
112     int solve(void* mem, const double* A, double* x, casadi_int nrhs, bool tr) const override;
113 
114     /// A documentation string
115     static const std::string meta_doc;
116 
117     // Get name of the plugin
plugin_name() const118     const char* plugin_name() const override { return "lapackqr";}
119 
120     // Get name of the class
class_name() const121     std::string class_name() const override { return "LapackQr";}
122 
123     // Maximum number of right-hand-sides
124     casadi_int max_nrhs_;
125 
126     /** \brief Serialize an object without type information */
127     void serialize_body(SerializingStream &s) const override;
128 
129     /** \brief Deserialize with type disambiguation */
deserialize(DeserializingStream & s)130     static ProtoFunction* deserialize(DeserializingStream& s) { return new LapackQr(s); }
131 
132   protected:
133     /** \brief Deserializing constructor */
134     explicit LapackQr(DeserializingStream& s);
135   };
136 
137 } // namespace casadi
138 
139 /// \endcond
140 #endif // CASADI_LAPACK_QR_HPP
141