dbleLU.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1994-2012 John W. Eaton
00004 Copyright (C) 2009 VZLU Prague, a.s.
00005 
00006 This file is part of Octave.
00007 
00008 Octave is free software; you can redistribute it and/or modify it
00009 under the terms of the GNU General Public License as published by the
00010 Free Software Foundation; either version 3 of the License, or (at your
00011 option) any later version.
00012 
00013 Octave is distributed in the hope that it will be useful, but WITHOUT
00014 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00015 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00016 for more details.
00017 
00018 You should have received a copy of the GNU General Public License
00019 along with Octave; see the file COPYING.  If not, see
00020 <http://www.gnu.org/licenses/>.
00021 
00022 */
00023 
00024 #ifdef HAVE_CONFIG_H
00025 #include <config.h>
00026 #endif
00027 
00028 #include "dbleLU.h"
00029 #include "f77-fcn.h"
00030 #include "lo-error.h"
00031 #include "oct-locbuf.h"
00032 #include "dColVector.h"
00033 
00034 // Instantiate the base LU class for the types we need.
00035 
00036 #include <base-lu.h>
00037 #include <base-lu.cc>
00038 
00039 template class base_lu <Matrix>;
00040 
00041 // Define the constructor for this particular derivation.
00042 
00043 extern "C"
00044 {
00045   F77_RET_T
00046   F77_FUNC (dgetrf, DGETRF) (const octave_idx_type&, const octave_idx_type&,
00047                              double*, const octave_idx_type&,
00048                              octave_idx_type*, octave_idx_type&);
00049 
00050 #ifdef HAVE_QRUPDATE_LUU
00051   F77_RET_T
00052   F77_FUNC (dlu1up, DLU1UP) (const octave_idx_type&, const octave_idx_type&,
00053                              double *, const octave_idx_type&,
00054                              double *, const octave_idx_type&,
00055                              double *, double *);
00056 
00057   F77_RET_T
00058   F77_FUNC (dlup1up, DLUP1UP) (const octave_idx_type&, const octave_idx_type&,
00059                                double *, const octave_idx_type&,
00060                                double *, const octave_idx_type&,
00061                                octave_idx_type *, const double *,
00062                                const double *, double *);
00063 #endif
00064 }
00065 
00066 LU::LU (const Matrix& a)
00067 {
00068   octave_idx_type a_nr = a.rows ();
00069   octave_idx_type a_nc = a.cols ();
00070   octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
00071 
00072   ipvt.resize (dim_vector (mn, 1));
00073   octave_idx_type *pipvt = ipvt.fortran_vec ();
00074 
00075   a_fact = a;
00076   double *tmp_data = a_fact.fortran_vec ();
00077 
00078   octave_idx_type info = 0;
00079 
00080   F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
00081 
00082   for (octave_idx_type i = 0; i < mn; i++)
00083     pipvt[i] -= 1;
00084 }
00085 
00086 #ifdef HAVE_QRUPDATE_LUU
00087 
00088 void LU::update (const ColumnVector& u, const ColumnVector& v)
00089 {
00090   if (packed ())
00091     unpack ();
00092 
00093   Matrix& l = l_fact;
00094   Matrix& r = a_fact;
00095 
00096   octave_idx_type m = l.rows ();
00097   octave_idx_type n = r.columns ();
00098   octave_idx_type k = l.columns ();
00099 
00100   if (u.length () == m && v.length () == n)
00101     {
00102       ColumnVector utmp = u, vtmp = v;
00103       F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
00104                                  utmp.fortran_vec (), vtmp.fortran_vec ()));
00105     }
00106   else
00107     (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
00108 }
00109 
00110 void LU::update (const Matrix& u, const Matrix& v)
00111 {
00112   if (packed ())
00113     unpack ();
00114 
00115   Matrix& l = l_fact;
00116   Matrix& r = a_fact;
00117 
00118   octave_idx_type m = l.rows ();
00119   octave_idx_type n = r.columns ();
00120   octave_idx_type k = l.columns ();
00121 
00122   if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
00123     {
00124       for (volatile octave_idx_type i = 0; i < u.cols (); i++)
00125         {
00126           ColumnVector utmp = u.column (i), vtmp = v.column (i);
00127           F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
00128                                      utmp.fortran_vec (), vtmp.fortran_vec ()));
00129         }
00130     }
00131   else
00132     (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
00133 }
00134 
00135 void LU::update_piv (const ColumnVector& u, const ColumnVector& v)
00136 {
00137   if (packed ())
00138     unpack ();
00139 
00140   Matrix& l = l_fact;
00141   Matrix& r = a_fact;
00142 
00143   octave_idx_type m = l.rows ();
00144   octave_idx_type n = r.columns ();
00145   octave_idx_type k = l.columns ();
00146 
00147   if (u.length () == m && v.length () == n)
00148     {
00149       ColumnVector utmp = u, vtmp = v;
00150       OCTAVE_LOCAL_BUFFER (double, w, m);
00151       for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
00152       F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
00153                                    ipvt.fortran_vec (), utmp.data (), vtmp.data (), w));
00154       for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
00155     }
00156   else
00157     (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
00158 }
00159 
00160 void LU::update_piv (const Matrix& u, const Matrix& v)
00161 {
00162   if (packed ())
00163     unpack ();
00164 
00165   Matrix& l = l_fact;
00166   Matrix& r = a_fact;
00167 
00168   octave_idx_type m = l.rows ();
00169   octave_idx_type n = r.columns ();
00170   octave_idx_type k = l.columns ();
00171 
00172   if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
00173     {
00174       OCTAVE_LOCAL_BUFFER (double, w, m);
00175       for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
00176       for (volatile octave_idx_type i = 0; i < u.cols (); i++)
00177         {
00178           ColumnVector utmp = u.column (i), vtmp = v.column (i);
00179           F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
00180                                        ipvt.fortran_vec (), utmp.data (), vtmp.data (), w));
00181         }
00182       for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
00183     }
00184   else
00185     (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
00186 }
00187 
00188 #else
00189 
00190 void LU::update (const ColumnVector&, const ColumnVector&)
00191 {
00192   (*current_liboctave_error_handler) ("luupdate: not available in this version");
00193 }
00194 
00195 void LU::update (const Matrix&, const Matrix&)
00196 {
00197   (*current_liboctave_error_handler) ("luupdate: not available in this version");
00198 }
00199 
00200 void LU::update_piv (const ColumnVector&, const ColumnVector&)
00201 {
00202   (*current_liboctave_error_handler) ("luupdate: not available in this version");
00203 }
00204 
00205 void LU::update_piv (const Matrix&, const Matrix&)
00206 {
00207   (*current_liboctave_error_handler) ("luupdate: not available in this version");
00208 }
00209 
00210 #endif
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines