Mercurial > hg > octave-nkf

comparison liboctave/CmplxLU.cc @ 1992:5668c00f9c7a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-03-03 00:40:53 by jwe]
author jwe
date Sun, 03 Mar 1996 00:43:53 +0000
parents 3ce2c289c978
children 1b57120c997b
comparison
equal deleted inserted replaced
1991:413f6a81868f 1992:5668c00f9c7a
32 #include "CmplxLU.h" 32 #include "CmplxLU.h"
33 #include "f77-fcn.h" 33 #include "f77-fcn.h"
34 #include "lo-error.h" 34 #include "lo-error.h"
35 #include "mx-inlines.cc" 35 #include "mx-inlines.cc"
36 36
37 // Instantiate the base LU class for the types we need.
38
39 #include <base-lu.h>
40 #include <base-lu.cc>
41
42 template class base_lu <ComplexMatrix, Matrix>;
43
44 // Define the constructor for this particular derivation.
45
37 extern "C" 46 extern "C"
38 { 47 {
39 int F77_FCN (zgesv, ZGESV) (const int&, const int&, Complex*, 48 int F77_FCN (zgesv, ZGESV) (const int&, const int&, Complex*,
40 const int&, int*, Complex*, const int&, 49 const int&, int*, Complex*, const int&,
41 int&); 50 int&);
52 return; 61 return;
53 } 62 }
54 63
55 int n = a_nr; 64 int n = a_nr;
56 65
57 Array<int> ipvt (n); 66 ipvt.resize (n);
58
59 int *pipvt = ipvt.fortran_vec (); 67 int *pipvt = ipvt.fortran_vec ();
60 68
61 ComplexMatrix A_fact = a; 69 a_fact = a;
62 Complex *tmp_data = A_fact.fortran_vec (); 70 Complex *tmp_data = a_fact.fortran_vec ();
63 71
64 int info = 0; 72 int info = 0;
65 Complex *dummy = 0; 73 Complex *dummy = 0;
66 74
67 F77_XFCN (zgesv, ZGESV, (n, 0, tmp_data, n, pipvt, dummy, n, info)); 75 F77_XFCN (zgesv, ZGESV, (n, 0, tmp_data, n, pipvt, dummy, n, info));
68 76
69 if (f77_exception_encountered) 77 if (f77_exception_encountered)
70 (*current_liboctave_error_handler) ("unrecoverable error in zgesv"); 78 (*current_liboctave_error_handler) ("unrecoverable error in zgesv");
71 else 79 else
72 { 80 ipvt -= 1;
73 Array<int> pvt (n);
74
75 for (int i = 0; i < n; i++)
76 {
77 ipvt.elem (i) -= 1;
78 pvt.elem (i) = i;
79 }
80
81 for (int i = 0; i < n - 1; i++)
82 {
83 int k = ipvt.elem (i);
84
85 if (k != i)
86 {
87 int tmp = pvt.elem (k);
88 pvt.elem (k) = pvt.elem (i);
89 pvt.elem (i)= tmp;
90 }
91 }
92
93 l.resize (n, n, 0.0);
94 u.resize (n, n, 0.0);
95 p.resize (n, n, 0.0);
96
97 for (int i = 0; i < n; i++)
98 {
99 p.elem (i, pvt.elem (i)) = 1.0;
100
101 l.elem (i, i) = 1.0;
102 for (int j = 0; j < i; j++)
103 l.elem (i, j) = A_fact.elem (i, j);
104
105 for (int j = i; j < n; j++)
106 u.elem (i, j) = A_fact.elem (i, j);
107 }
108 }
109 } 81 }
110 82
111 /* 83 /*
112 ;;; Local Variables: *** 84 ;;; Local Variables: ***
113 ;;; mode: C++ *** 85 ;;; mode: C++ ***