00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00025 #pragma implementation
00026 #endif
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032 #include <iostream>
00033
00034 #include "Array-util.h"
00035 #include "f77-fcn.h"
00036 #include "lo-error.h"
00037 #include "mx-base.h"
00038 #include "mx-inlines.cc"
00039 #include "oct-cmplx.h"
00040
00041
00042
00043 extern "C"
00044 {
00045 F77_RET_T
00046 F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL,
00047 const int&, const int&, const double&,
00048 const double*, const int&, const double*,
00049 const int&, const double&, double*,
00050 const int&
00051 F77_CHAR_ARG_LEN_DECL);
00052 }
00053
00054
00055
00056 bool
00057 ColumnVector::operator == (const ColumnVector& a) const
00058 {
00059 int len = length ();
00060 if (len != a.length ())
00061 return 0;
00062 return mx_inline_equal (data (), a.data (), len);
00063 }
00064
00065 bool
00066 ColumnVector::operator != (const ColumnVector& a) const
00067 {
00068 return !(*this == a);
00069 }
00070
00071 ColumnVector&
00072 ColumnVector::insert (const ColumnVector& a, int r)
00073 {
00074 int a_len = a.length ();
00075
00076 if (r < 0 || r + a_len > length ())
00077 {
00078 (*current_liboctave_error_handler) ("range error for insert");
00079 return *this;
00080 }
00081
00082 if (a_len > 0)
00083 {
00084 make_unique ();
00085
00086 for (int i = 0; i < a_len; i++)
00087 xelem (r+i) = a.elem (i);
00088 }
00089
00090 return *this;
00091 }
00092
00093 ColumnVector&
00094 ColumnVector::fill (double val)
00095 {
00096 int len = length ();
00097
00098 if (len > 0)
00099 {
00100 make_unique ();
00101
00102 for (int i = 0; i < len; i++)
00103 xelem (i) = val;
00104 }
00105
00106 return *this;
00107 }
00108
00109 ColumnVector&
00110 ColumnVector::fill (double val, int r1, int r2)
00111 {
00112 int len = length ();
00113
00114 if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len)
00115 {
00116 (*current_liboctave_error_handler) ("range error for fill");
00117 return *this;
00118 }
00119
00120 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00121
00122 if (r2 >= r1)
00123 {
00124 make_unique ();
00125
00126 for (int i = r1; i <= r2; i++)
00127 xelem (i) = val;
00128 }
00129
00130 return *this;
00131 }
00132
00133 ColumnVector
00134 ColumnVector::stack (const ColumnVector& a) const
00135 {
00136 int len = length ();
00137 int nr_insert = len;
00138 ColumnVector retval (len + a.length ());
00139 retval.insert (*this, 0);
00140 retval.insert (a, nr_insert);
00141 return retval;
00142 }
00143
00144 RowVector
00145 ColumnVector::transpose (void) const
00146 {
00147 return RowVector (*this);
00148 }
00149
00150 ColumnVector
00151 real (const ComplexColumnVector& a)
00152 {
00153 int a_len = a.length ();
00154 ColumnVector retval;
00155 if (a_len > 0)
00156 retval = ColumnVector (mx_inline_real_dup (a.data (), a_len), a_len);
00157 return retval;
00158 }
00159
00160 ColumnVector
00161 imag (const ComplexColumnVector& a)
00162 {
00163 int a_len = a.length ();
00164 ColumnVector retval;
00165 if (a_len > 0)
00166 retval = ColumnVector (mx_inline_imag_dup (a.data (), a_len), a_len);
00167 return retval;
00168 }
00169
00170
00171
00172 ColumnVector
00173 ColumnVector::extract (int r1, int r2) const
00174 {
00175 if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
00176
00177 int new_r = r2 - r1 + 1;
00178
00179 ColumnVector result (new_r);
00180
00181 for (int i = 0; i < new_r; i++)
00182 result.xelem (i) = elem (r1+i);
00183
00184 return result;
00185 }
00186
00187 ColumnVector
00188 ColumnVector::extract_n (int r1, int n) const
00189 {
00190 ColumnVector result (n);
00191
00192 for (int i = 0; i < n; i++)
00193 result.xelem (i) = elem (r1+i);
00194
00195 return result;
00196 }
00197
00198
00199
00200 ColumnVector
00201 operator * (const Matrix& m, const ColumnVector& a)
00202 {
00203 ColumnVector retval;
00204
00205 int nr = m.rows ();
00206 int nc = m.cols ();
00207
00208 int a_len = a.length ();
00209
00210 if (nc != a_len)
00211 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
00212 else
00213 {
00214 if (nr == 0 || nc == 0)
00215 retval.resize (nr, 0.0);
00216 else
00217 {
00218 int ld = nr;
00219
00220 retval.resize (nr);
00221 double *y = retval.fortran_vec ();
00222
00223 F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
00224 nr, nc, 1.0, m.data (), ld,
00225 a.data (), 1, 0.0, y, 1
00226 F77_CHAR_ARG_LEN (1)));
00227
00228 if (f77_exception_encountered)
00229 (*current_liboctave_error_handler)
00230 ("unrecoverable error in dgemv");
00231 }
00232 }
00233
00234 return retval;
00235 }
00236
00237
00238
00239 ColumnVector
00240 operator * (const DiagMatrix& m, const ColumnVector& a)
00241 {
00242 ColumnVector retval;
00243
00244 int nr = m.rows ();
00245 int nc = m.cols ();
00246
00247 int a_len = a.length ();
00248
00249 if (nc != a_len)
00250 gripe_nonconformant ("operator *", nr, nc, a_len, 1);
00251 else
00252 {
00253 if (nr == 0 || nc == 0)
00254 retval.resize (nr, 0.0);
00255 else
00256 {
00257 retval.resize (nr);
00258
00259 for (int i = 0; i < a_len; i++)
00260 retval.elem (i) = a.elem (i) * m.elem (i, i);
00261
00262 for (int i = a_len; i < nr; i++)
00263 retval.elem (i) = 0.0;
00264 }
00265 }
00266
00267 return retval;
00268 }
00269
00270
00271
00272 ColumnVector
00273 ColumnVector::map (d_d_Mapper f) const
00274 {
00275 ColumnVector b (*this);
00276 return b.apply (f);
00277 }
00278
00279 ColumnVector&
00280 ColumnVector::apply (d_d_Mapper f)
00281 {
00282 double *d = fortran_vec ();
00283
00284 for (int i = 0; i < length (); i++)
00285 d[i] = f (d[i]);
00286
00287 return *this;
00288 }
00289
00290 double
00291 ColumnVector::min (void) const
00292 {
00293 int len = length ();
00294 if (len == 0)
00295 return 0.0;
00296
00297 double res = elem (0);
00298
00299 for (int i = 1; i < len; i++)
00300 if (elem (i) < res)
00301 res = elem (i);
00302
00303 return res;
00304 }
00305
00306 double
00307 ColumnVector::max (void) const
00308 {
00309 int len = length ();
00310 if (len == 0)
00311 return 0.0;
00312
00313 double res = elem (0);
00314
00315 for (int i = 1; i < len; i++)
00316 if (elem (i) > res)
00317 res = elem (i);
00318
00319 return res;
00320 }
00321
00322 std::ostream&
00323 operator << (std::ostream& os, const ColumnVector& a)
00324 {
00325
00326 for (int i = 0; i < a.length (); i++)
00327 os << a.elem (i) << "\n";
00328 return os;
00329 }
00330
00331 std::istream&
00332 operator >> (std::istream& is, ColumnVector& a)
00333 {
00334 int len = a.length();
00335
00336 if (len < 1)
00337 is.clear (std::ios::badbit);
00338 else
00339 {
00340 double tmp;
00341 for (int i = 0; i < len; i++)
00342 {
00343 is >> tmp;
00344 if (is)
00345 a.elem (i) = tmp;
00346 else
00347 break;
00348 }
00349 }
00350 return is;
00351 }
00352
00353
00354
00355
00356
00357