00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00024 #pragma implementation
00025 #endif
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif
00030
00031 #include "dbleCHOL.h"
00032 #include "f77-fcn.h"
00033 #include "lo-error.h"
00034
00035 extern "C"
00036 {
00037 F77_RET_T
00038 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, const int&,
00039 double*, const int&, int&
00040 F77_CHAR_ARG_LEN_DECL);
00041 }
00042
00043 int
00044 CHOL::init (const Matrix& a)
00045 {
00046 int a_nr = a.rows ();
00047 int a_nc = a.cols ();
00048
00049 if (a_nr != a_nc)
00050 {
00051 (*current_liboctave_error_handler) ("CHOL requires square matrix");
00052 return -1;
00053 }
00054
00055 int n = a_nc;
00056 int info;
00057
00058 chol_mat = a;
00059 double *h = chol_mat.fortran_vec ();
00060
00061 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
00062 n, h, n, info
00063 F77_CHAR_ARG_LEN (1)));
00064
00065 if (f77_exception_encountered)
00066 (*current_liboctave_error_handler) ("unrecoverable error in dpotrf");
00067 else
00068 {
00069
00070
00071
00072 if (n > 1)
00073 for (int j = 0; j < a_nc; j++)
00074 for (int i = j+1; i < a_nr; i++)
00075 chol_mat.elem (i, j) = 0.0;
00076 }
00077
00078 return info;
00079 }
00080
00081
00082
00083
00084
00085