Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/matrix_type.cc @ 5785:6b9cec830d72

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-05-03 19:32:46 by dbateman]
author dbateman
date Wed, 03 May 2006 19:32:48 +0000
parents ace8d8d26933
children 080c08b192d8
comparison
equal deleted inserted replaced
5784:70f7659d0fb9 5785:6b9cec830d72
28 #include <algorithm> 28 #include <algorithm>
29 29
30 #include "ov.h" 30 #include "ov.h"
31 #include "defun-dld.h" 31 #include "defun-dld.h"
32 #include "error.h" 32 #include "error.h"
33 #include "ov-re-mat.h"
34 #include "ov-cx-mat.h"
33 #include "ov-re-sparse.h" 35 #include "ov-re-sparse.h"
34 #include "ov-cx-sparse.h" 36 #include "ov-cx-sparse.h"
35 #include "SparseType.h" 37 #include "MatrixType.h"
36 38
37 DEFUN_DLD (matrix_type, args, , 39 DEFUN_DLD (matrix_type, args, ,
38 "-*- texinfo -*-\n\ 40 "-*- texinfo -*-\n\
39 @deftypefn {Loadable Function} {@var{type} =} matrix_type (@var{a})\n\ 41 @deftypefn {Loadable Function} {@var{type} =} matrix_type (@var{a})\n\
40 @deftypefnx {Loadable Function} {@var{a} =} matrix_type (@var{a}, @var{type})\n\ 42 @deftypefnx {Loadable Function} {@var{a} =} matrix_type (@var{a}, @var{type})\n\
109 { 111 {
110 if (args(0).is_sparse_type ()) 112 if (args(0).is_sparse_type ())
111 { 113 {
112 if (nargin == 1) 114 if (nargin == 1)
113 { 115 {
114 SparseType mattyp; 116 MatrixType mattyp;
115 117
116 if (args(0).type_name () == "sparse complex matrix" ) 118 if (args(0).type_name () == "sparse complex matrix" )
117 { 119 {
118 const octave_sparse_complex_matrix& rep 120 mattyp = args(0).matrix_type ();
119 = dynamic_cast<const octave_sparse_complex_matrix&> (args(0).get_rep ());
120
121 mattyp = rep.sparse_type ();
122 121
123 if (mattyp.is_unknown ()) 122 if (mattyp.is_unknown ())
124 { 123 {
125 SparseComplexMatrix m = 124 SparseComplexMatrix m =
126 args(0).sparse_complex_matrix_value (); 125 args(0).sparse_complex_matrix_value ();
127 if (!error_state) 126 if (!error_state)
128 { 127 {
129 mattyp = SparseType (m); 128 mattyp = MatrixType (m);
130 rep.sparse_type (mattyp); 129 args(0).matrix_type (mattyp);
131 } 130 }
132 } 131 }
133 } 132 }
134 else 133 else
135 { 134 {
136 const octave_sparse_matrix& rep 135 mattyp = args(0).matrix_type ();
137 = dynamic_cast<const octave_sparse_matrix&> (args(0).get_rep ());
138
139 mattyp = rep.sparse_type ();
140 136
141 if (mattyp.is_unknown ()) 137 if (mattyp.is_unknown ())
142 { 138 {
143 SparseMatrix m = args(0).sparse_matrix_value (); 139 SparseMatrix m = args(0).sparse_matrix_value ();
144 if (!error_state) 140 if (!error_state)
145 { 141 {
146 mattyp = SparseType (m); 142 mattyp = MatrixType (m);
147 rep.sparse_type (mattyp); 143 args(0).matrix_type (mattyp);
148 } 144 }
149 } 145 }
150 } 146 }
151 147
152 int typ = mattyp.type (); 148 int typ = mattyp.type ();
153 149
154 if (typ == SparseType::Diagonal) 150 if (typ == MatrixType::Diagonal)
155 retval = octave_value ("Diagonal"); 151 retval = octave_value ("Diagonal");
156 else if (typ == SparseType::Permuted_Diagonal) 152 else if (typ == MatrixType::Permuted_Diagonal)
157 retval = octave_value ("Permuted Diagonal"); 153 retval = octave_value ("Permuted Diagonal");
158 else if (typ == SparseType::Upper) 154 else if (typ == MatrixType::Upper)
159 retval = octave_value ("Upper"); 155 retval = octave_value ("Upper");
160 else if (typ == SparseType::Permuted_Upper) 156 else if (typ == MatrixType::Permuted_Upper)
161 retval = octave_value ("Permuted Upper"); 157 retval = octave_value ("Permuted Upper");
162 else if (typ == SparseType::Lower) 158 else if (typ == MatrixType::Lower)
163 retval = octave_value ("Lower"); 159 retval = octave_value ("Lower");
164 else if (typ == SparseType::Permuted_Lower) 160 else if (typ == MatrixType::Permuted_Lower)
165 retval = octave_value ("Permuted Lower"); 161 retval = octave_value ("Permuted Lower");
166 else if (typ == SparseType::Banded) 162 else if (typ == MatrixType::Banded)
167 retval = octave_value ("Banded"); 163 retval = octave_value ("Banded");
168 else if (typ == SparseType::Banded_Hermitian) 164 else if (typ == MatrixType::Banded_Hermitian)
169 retval = octave_value ("Banded Positive Definite"); 165 retval = octave_value ("Banded Positive Definite");
170 else if (typ == SparseType::Tridiagonal) 166 else if (typ == MatrixType::Tridiagonal)
171 retval = octave_value ("Tridiagonal"); 167 retval = octave_value ("Tridiagonal");
172 else if (typ == SparseType::Tridiagonal_Hermitian) 168 else if (typ == MatrixType::Tridiagonal_Hermitian)
173 retval = octave_value ("Tridiagonal Positive Definite"); 169 retval = octave_value ("Tridiagonal Positive Definite");
174 else if (typ == SparseType::Hermitian) 170 else if (typ == MatrixType::Hermitian)
175 retval = octave_value ("Positive Definite"); 171 retval = octave_value ("Positive Definite");
176 else if (typ == SparseType::Rectangular) 172 else if (typ == MatrixType::Rectangular)
177 { 173 {
178 if (args(0).rows() == args(0).columns()) 174 if (args(0).rows() == args(0).columns())
179 retval = octave_value ("Singular"); 175 retval = octave_value ("Singular");
180 else 176 else
181 retval = octave_value ("Rectangular"); 177 retval = octave_value ("Rectangular");
182 } 178 }
183 else if (typ == SparseType::Full) 179 else if (typ == MatrixType::Full)
184 retval = octave_value ("Full"); 180 retval = octave_value ("Full");
185 else 181 else
186 // This should never happen!!! 182 // This should never happen!!!
187 retval = octave_value ("Unknown"); 183 retval = octave_value ("Unknown");
188 } 184 }
190 { 186 {
191 // Ok, we're changing the matrix type 187 // Ok, we're changing the matrix type
192 std::string str_typ = args(1).string_value (); 188 std::string str_typ = args(1).string_value ();
193 189
194 // FIXME -- why do I have to explicitly call the constructor? 190 // FIXME -- why do I have to explicitly call the constructor?
195 SparseType mattyp = SparseType (); 191 MatrixType mattyp = MatrixType ();
196 192
197 octave_idx_type nl = 0; 193 octave_idx_type nl = 0;
198 octave_idx_type nu = 0; 194 octave_idx_type nu = 0;
199 195
200 if (error_state) 196 if (error_state)
298 } 294 }
299 } 295 }
300 } 296 }
301 } 297 }
302 else 298 else
303 error ("matrix_type: Only sparse matrices treated at the moment"); 299 {
300 if (nargin == 1)
301 {
302 MatrixType mattyp;
303
304 if (args(0).type_name () == "complex matrix" )
305 {
306 mattyp = args(0).matrix_type ();
307
308 if (mattyp.is_unknown ())
309 {
310 ComplexMatrix m = args(0).complex_matrix_value ();
311 if (!error_state)
312 {
313 mattyp = MatrixType (m);
314 args(0).matrix_type (mattyp);
315 }
316 }
317 }
318 else
319 {
320 mattyp = args(0).matrix_type ();
321
322 if (mattyp.is_unknown ())
323 {
324 Matrix m = args(0).matrix_value ();
325 if (!error_state)
326 {
327 mattyp = MatrixType (m);
328 args(0).matrix_type (mattyp);
329 }
330 }
331 }
332
333 int typ = mattyp.type ();
334
335 if (typ == MatrixType::Upper)
336 retval = octave_value ("Upper");
337 else if (typ == MatrixType::Permuted_Upper)
338 retval = octave_value ("Permuted Upper");
339 else if (typ == MatrixType::Lower)
340 retval = octave_value ("Lower");
341 else if (typ == MatrixType::Permuted_Lower)
342 retval = octave_value ("Permuted Lower");
343 else if (typ == MatrixType::Hermitian)
344 retval = octave_value ("Positive Definite");
345 else if (typ == MatrixType::Rectangular)
346 {
347 if (args(0).rows() == args(0).columns())
348 retval = octave_value ("Singular");
349 else
350 retval = octave_value ("Rectangular");
351 }
352 else if (typ == MatrixType::Full)
353 retval = octave_value ("Full");
354 else
355 // This should never happen!!!
356 retval = octave_value ("Unknown");
357 }
358 else
359 {
360 // Ok, we're changing the matrix type
361 std::string str_typ = args(1).string_value ();
362
363 // FIXME -- why do I have to explicitly call the constructor?
364 MatrixType mattyp = MatrixType (MatrixType::Unknown, true);
365
366 if (error_state)
367 error ("Matrix type must be a string");
368 else
369 {
370 // Use STL function to convert to lower case
371 std::transform (str_typ.begin (), str_typ.end (),
372 str_typ.begin (), tolower);
373
374 if (str_typ == "upper")
375 mattyp.mark_as_upper_triangular ();
376 else if (str_typ == "lower")
377 mattyp.mark_as_lower_triangular ();
378 else if (str_typ == "positive definite")
379 {
380 mattyp.mark_as_full ();
381 mattyp.mark_as_symmetric ();
382 }
383 else if (str_typ == "singular")
384 mattyp.mark_as_rectangular ();
385 else if (str_typ == "full")
386 mattyp.mark_as_full ();
387 else if (str_typ == "unknown")
388 mattyp.invalidate_type ();
389 else
390 error ("matrix_type: Unknown matrix type %s", str_typ.c_str());
391
392 if (! error_state)
393 {
394 if (nargin == 3 && (str_typ == "upper"
395 || str_typ == "lower"))
396 {
397 const ColumnVector perm =
398 ColumnVector (args (2).vector_value ());
399
400 if (error_state)
401 error ("matrix_type: Invalid permutation vector");
402 else
403 {
404 octave_idx_type len = perm.length ();
405 dim_vector dv = args(0).dims ();
406
407 if (len != dv(0))
408 error ("matrix_type: Invalid permutation vector");
409 else
410 {
411 OCTAVE_LOCAL_BUFFER (octave_idx_type, p, len);
412
413 for (octave_idx_type i = 0; i < len; i++)
414 p[i] = static_cast<octave_idx_type> (perm (i)) - 1;
415
416 if (str_typ == "upper")
417 mattyp.mark_as_permuted (len, p);
418 else
419 mattyp.mark_as_permuted (len, p);
420 }
421 }
422 }
423 else if (nargin != 2)
424 error ("matrix_type: Invalid number of arguments");
425
426 if (! error_state)
427 {
428 // Set the matrix type
429 if (args(0).type_name () == "complex matrix" )
430 retval =
431 octave_value (args(0).complex_matrix_value (),
432 mattyp);
433 else
434 retval = octave_value (args(0).matrix_value (),
435 mattyp);
436 }
437 }
438 }
439 }
440 }
304 } 441 }
305 442
306 return retval; 443 return retval;
307 } 444 }
308 445
309 /* 446 /*
310 447
311 ## FIXME 448 ## FIXME
312 ## Disable tests for lower under-determined and upper over-determined 449 ## Disable tests for lower under-determined and upper over-determined
313 ## matrices and this detection is disabled in SparseType due to issues 450 ## matrices and this detection is disabled in MatrixType due to issues
314 ## of non minimum norm solution being found. 451 ## of non minimum norm solution being found.
315 452
316 %!assert(matrix_type(speye(10,10)),"Diagonal"); 453 %!assert(matrix_type(speye(10,10)),"Diagonal");
317 %!assert(matrix_type(speye(10,10)([2:10,1],:)),"Permuted Diagonal"); 454 %!assert(matrix_type(speye(10,10)([2:10,1],:)),"Permuted Diagonal");
318 %!assert(matrix_type([[speye(10,10);sparse(1,10)],[1;sparse(9,1);1]]),"Upper"); 455 %!assert(matrix_type([[speye(10,10);sparse(1,10)],[1;sparse(9,1);1]]),"Upper");
391 528
392 %!test 529 %!test
393 %! a = matrix_type(spdiags(randn(10,3),[-1,0,1],10,10),"Singular"); 530 %! a = matrix_type(spdiags(randn(10,3),[-1,0,1],10,10),"Singular");
394 %! assert(matrix_type(a),"Singular"); 531 %! assert(matrix_type(a),"Singular");
395 532
533 %!assert(matrix_type(triu(ones(10,10))),"Upper");
534 %!assert(matrix_type(triu(ones(10,10),-1)),"Full");
535 %!assert(matrix_type(tril(ones(10,10))),"Lower");
536 %!assert(matrix_type(tril(ones(10,10),1)),"Full");
537 %!assert(matrix_type(10*eye(10,10) + ones(10,10)), "Positive Definite");
538 %!assert(matrix_type(ones(11,10)),"Rectangular")
539 %!test
540 %! a = matrix_type(ones(10,10),"Singular");
541 %! assert(matrix_type(a),"Singular");
542
543 %!assert(matrix_type(triu(1i*ones(10,10))),"Upper");
544 %!assert(matrix_type(triu(1i*ones(10,10),-1)),"Full");
545 %!assert(matrix_type(tril(1i*ones(10,10))),"Lower");
546 %!assert(matrix_type(tril(1i*ones(10,10),1)),"Full");
547 %!assert(matrix_type(10*eye(10,10) + 1i*triu(ones(10,10),1) -1i*tril(ones(10,10),-1)), "Positive Definite");
548 %!assert(matrix_type(ones(11,10)),"Rectangular")
549 %!test
550 %! a = matrix_type(ones(10,10),"Singular");
551 %! assert(matrix_type(a),"Singular");
552
396 */ 553 */
397 554
398 /* 555 /*
399 ;;; Local Variables: *** 556 ;;; Local Variables: ***
400 ;;; mode: C++ *** 557 ;;; mode: C++ ***