Mercurial > hg > octave-nkf

comparison liboctave/CMatrix.cc @ 9661:afcf852256d2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
optimize / and '\ for triangular matrices
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 23 Sep 2009 10:00:16 +0200
parents 3429c956de6f
children 7e5b4de5fbfe
comparison
equal deleted inserted replaced
9660:0256e187d13b 9661:afcf852256d2
1850 1850
1851 ComplexMatrix 1851 ComplexMatrix
1852 ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b, 1852 ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
1853 octave_idx_type& info, double& rcon, 1853 octave_idx_type& info, double& rcon,
1854 solve_singularity_handler sing_handler, 1854 solve_singularity_handler sing_handler,
1855 bool calc_cond) const 1855 bool calc_cond, blas_trans_type transt) const
1856 { 1856 {
1857 ComplexMatrix retval; 1857 ComplexMatrix retval;
1858 1858
1859 octave_idx_type nr = rows (); 1859 octave_idx_type nr = rows ();
1860 octave_idx_type nc = cols (); 1860 octave_idx_type nc = cols ();
1926 { 1926 {
1927 retval = b; 1927 retval = b;
1928 Complex *result = retval.fortran_vec (); 1928 Complex *result = retval.fortran_vec ();
1929 1929
1930 char uplo = 'U'; 1930 char uplo = 'U';
1931 char trans = 'N'; 1931 char trans = get_blas_char (transt);
1932 char dia = 'N'; 1932 char dia = 'N';
1933 1933
1934 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), 1934 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1935 F77_CONST_CHAR_ARG2 (&trans, 1), 1935 F77_CONST_CHAR_ARG2 (&trans, 1),
1936 F77_CONST_CHAR_ARG2 (&dia, 1), 1936 F77_CONST_CHAR_ARG2 (&dia, 1),
1951 1951
1952 ComplexMatrix 1952 ComplexMatrix
1953 ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b, 1953 ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
1954 octave_idx_type& info, double& rcon, 1954 octave_idx_type& info, double& rcon,
1955 solve_singularity_handler sing_handler, 1955 solve_singularity_handler sing_handler,
1956 bool calc_cond) const 1956 bool calc_cond, blas_trans_type transt) const
1957 { 1957 {
1958 ComplexMatrix retval; 1958 ComplexMatrix retval;
1959 1959
1960 octave_idx_type nr = rows (); 1960 octave_idx_type nr = rows ();
1961 octave_idx_type nc = cols (); 1961 octave_idx_type nc = cols ();
2027 { 2027 {
2028 retval = b; 2028 retval = b;
2029 Complex *result = retval.fortran_vec (); 2029 Complex *result = retval.fortran_vec ();
2030 2030
2031 char uplo = 'L'; 2031 char uplo = 'L';
2032 char trans = 'N'; 2032 char trans = get_blas_char (transt);
2033 char dia = 'N'; 2033 char dia = 'N';
2034 2034
2035 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), 2035 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
2036 F77_CONST_CHAR_ARG2 (&trans, 1), 2036 F77_CONST_CHAR_ARG2 (&trans, 1),
2037 F77_CONST_CHAR_ARG2 (&dia, 1), 2037 F77_CONST_CHAR_ARG2 (&dia, 1),
2258 } 2258 }
2259 2259
2260 ComplexMatrix 2260 ComplexMatrix
2261 ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, 2261 ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
2262 double& rcon, solve_singularity_handler sing_handler, 2262 double& rcon, solve_singularity_handler sing_handler,
2263 bool singular_fallback) const 2263 bool singular_fallback, blas_trans_type transt) const
2264 { 2264 {
2265 ComplexMatrix tmp (b); 2265 ComplexMatrix tmp (b);
2266 return solve (typ, tmp, info, rcon, sing_handler, singular_fallback); 2266 return solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt);
2267 } 2267 }
2268 2268
2269 ComplexMatrix 2269 ComplexMatrix
2270 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b) const 2270 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b) const
2271 { 2271 {
2291 2291
2292 ComplexMatrix 2292 ComplexMatrix
2293 ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b, 2293 ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
2294 octave_idx_type& info, double& rcon, 2294 octave_idx_type& info, double& rcon,
2295 solve_singularity_handler sing_handler, 2295 solve_singularity_handler sing_handler,
2296 bool singular_fallback) const 2296 bool singular_fallback, blas_trans_type transt) const
2297 { 2297 {
2298 ComplexMatrix retval; 2298 ComplexMatrix retval;
2299 int typ = mattype.type (); 2299 int typ = mattype.type ();
2300 2300
2301 if (typ == MatrixType::Unknown) 2301 if (typ == MatrixType::Unknown)
2302 typ = mattype.type (*this); 2302 typ = mattype.type (*this);
2303 2303
2304 // Only calculate the condition number for LU/Cholesky 2304 // Only calculate the condition number for LU/Cholesky
2305 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) 2305 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
2306 retval = utsolve (mattype, b, info, rcon, sing_handler, false); 2306 retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt);
2307 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) 2307 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
2308 retval = ltsolve (mattype, b, info, rcon, sing_handler, false); 2308 retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt);
2309 else if (transt == blas_trans)
2310 return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback);
2311 else if (transt == blas_conj_trans)
2312 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler, singular_fallback);
2309 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 2313 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
2310 retval = fsolve (mattype, b, info, rcon, sing_handler, true); 2314 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
2311 else if (typ != MatrixType::Rectangular) 2315 else if (typ != MatrixType::Rectangular)
2312 { 2316 {
2313 (*current_liboctave_error_handler) ("unknown matrix type"); 2317 (*current_liboctave_error_handler) ("unknown matrix type");
2348 } 2352 }
2349 2353
2350 ComplexColumnVector 2354 ComplexColumnVector
2351 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, 2355 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
2352 octave_idx_type& info, double& rcon, 2356 octave_idx_type& info, double& rcon,
2353 solve_singularity_handler sing_handler) const 2357 solve_singularity_handler sing_handler, blas_trans_type transt) const
2354 { 2358 {
2355 return solve (typ, ComplexColumnVector (b), info, rcon, sing_handler); 2359 return solve (typ, ComplexColumnVector (b), info, rcon, sing_handler, transt);
2356 } 2360 }
2357 2361
2358 ComplexColumnVector 2362 ComplexColumnVector
2359 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b) const 2363 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b) const
2360 { 2364 {
2379 } 2383 }
2380 2384
2381 ComplexColumnVector 2385 ComplexColumnVector
2382 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, 2386 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
2383 octave_idx_type& info, double& rcon, 2387 octave_idx_type& info, double& rcon,
2384 solve_singularity_handler sing_handler) const 2388 solve_singularity_handler sing_handler, blas_trans_type transt) const
2385 { 2389 {
2386 2390
2387 ComplexMatrix tmp (b); 2391 ComplexMatrix tmp (b);
2388 return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); 2392 return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0));
2389 } 2393 }
2390 2394
2391 ComplexMatrix 2395 ComplexMatrix
2392 ComplexMatrix::solve (const Matrix& b) const 2396 ComplexMatrix::solve (const Matrix& b) const
2393 { 2397 {
2409 return solve (b, info, rcon, 0); 2413 return solve (b, info, rcon, 0);
2410 } 2414 }
2411 2415
2412 ComplexMatrix 2416 ComplexMatrix
2413 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon, 2417 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon,
2414 solve_singularity_handler sing_handler) const 2418 solve_singularity_handler sing_handler, blas_trans_type transt) const
2415 { 2419 {
2416 ComplexMatrix tmp (b); 2420 ComplexMatrix tmp (b);
2417 return solve (tmp, info, rcon, sing_handler); 2421 return solve (tmp, info, rcon, sing_handler, transt);
2418 } 2422 }
2419 2423
2420 ComplexMatrix 2424 ComplexMatrix
2421 ComplexMatrix::solve (const ComplexMatrix& b) const 2425 ComplexMatrix::solve (const ComplexMatrix& b) const
2422 { 2426 {
2438 return solve (b, info, rcon, 0); 2442 return solve (b, info, rcon, 0);
2439 } 2443 }
2440 2444
2441 ComplexMatrix 2445 ComplexMatrix
2442 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, 2446 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon,
2443 solve_singularity_handler sing_handler) const 2447 solve_singularity_handler sing_handler, blas_trans_type transt) const
2444 { 2448 {
2445 MatrixType mattype (*this); 2449 MatrixType mattype (*this);
2446 return solve (mattype, b, info, rcon, sing_handler); 2450 return solve (mattype, b, info, rcon, sing_handler, true, transt);
2447 } 2451 }
2448 2452
2449 ComplexColumnVector 2453 ComplexColumnVector
2450 ComplexMatrix::solve (const ColumnVector& b) const 2454 ComplexMatrix::solve (const ColumnVector& b) const
2451 { 2455 {
2469 } 2473 }
2470 2474
2471 ComplexColumnVector 2475 ComplexColumnVector
2472 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, 2476 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
2473 double& rcon, 2477 double& rcon,
2474 solve_singularity_handler sing_handler) const 2478 solve_singularity_handler sing_handler, blas_trans_type transt) const
2475 { 2479 {
2476 return solve (ComplexColumnVector (b), info, rcon, sing_handler); 2480 return solve (ComplexColumnVector (b), info, rcon, sing_handler, transt);
2477 } 2481 }
2478 2482
2479 ComplexColumnVector 2483 ComplexColumnVector
2480 ComplexMatrix::solve (const ComplexColumnVector& b) const 2484 ComplexMatrix::solve (const ComplexColumnVector& b) const
2481 { 2485 {
2499 } 2503 }
2500 2504
2501 ComplexColumnVector 2505 ComplexColumnVector
2502 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, 2506 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
2503 double& rcon, 2507 double& rcon,
2504 solve_singularity_handler sing_handler) const 2508 solve_singularity_handler sing_handler, blas_trans_type transt) const
2505 { 2509 {
2506 MatrixType mattype (*this); 2510 MatrixType mattype (*this);
2507 return solve (mattype, b, info, rcon, sing_handler); 2511 return solve (mattype, b, info, rcon, sing_handler, transt);
2508 } 2512 }
2509 2513
2510 ComplexMatrix 2514 ComplexMatrix
2511 ComplexMatrix::lssolve (const Matrix& b) const 2515 ComplexMatrix::lssolve (const Matrix& b) const
2512 { 2516 {