Mercurial > hg > octave-nkf

comparison liboctave/CMatrix.cc @ 7071:c3b479e753dd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-26 15:14:34 by jwe]
author jwe
date Fri, 26 Oct 2007 15:14:35 +0000
parents f0142f2afdc6
children b48d486f641d
comparison
equal deleted inserted replaced
7070:7593f8e83a2e 7071:c3b479e753dd
118 const octave_idx_type&, const double&, double&, 118 const octave_idx_type&, const double&, double&,
119 Complex*, double*, octave_idx_type& 119 Complex*, double*, octave_idx_type&
120 F77_CHAR_ARG_LEN_DECL); 120 F77_CHAR_ARG_LEN_DECL);
121 121
122 F77_RET_T 122 F77_RET_T
123 F77_FUNC (zgelsy, ZGELSY) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 123 F77_FUNC (zgelss, ZGELSS) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
124 Complex*, const octave_idx_type&, Complex*, 124 Complex*, const octave_idx_type&, Complex*,
125 const octave_idx_type&, octave_idx_type*, double&, octave_idx_type&, 125 const octave_idx_type&, double*, double&, octave_idx_type&,
126 Complex*, const octave_idx_type&, double*, octave_idx_type&); 126 Complex*, const octave_idx_type&, double*, octave_idx_type&);
127 127
128 F77_RET_T 128 F77_RET_T
129 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, 129 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
130 Complex*, const octave_idx_type&, 130 Complex*, const octave_idx_type&,
2446 for (octave_idx_type i = 0; i < m; i++) 2446 for (octave_idx_type i = 0; i < m; i++)
2447 result.elem (i, j) = b.elem (i, j); 2447 result.elem (i, j) = b.elem (i, j);
2448 2448
2449 Complex *presult = result.fortran_vec (); 2449 Complex *presult = result.fortran_vec ();
2450 2450
2451 Array<octave_idx_type> jpvt (n); 2451 octave_idx_type len_s = m < n ? m : n;
2452 octave_idx_type *pjpvt = jpvt.fortran_vec (); 2452 Array<double> s (len_s);
2453 double *ps = s.fortran_vec ();
2453 2454
2454 double rcond = -1.0; 2455 double rcond = -1.0;
2455 2456
2456 Array<double> rwork (2 * n); 2457 octave_idx_type lrwork = (5 * (m < n ? m : n)) - 4;
2458 lrwork = lrwork > 1 ? lrwork : 1;
2459 Array<double> rwork (lrwork);
2457 double *prwork = rwork.fortran_vec (); 2460 double *prwork = rwork.fortran_vec ();
2458 2461
2459 // Ask ZGELSY what the dimension of WORK should be. 2462 // Ask ZGELSS what the dimension of WORK should be.
2460 2463
2461 octave_idx_type lwork = -1; 2464 octave_idx_type lwork = -1;
2462 2465
2463 Array<Complex> work (1); 2466 Array<Complex> work (1);
2464 2467
2465 F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, 2468 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
2466 nrr, pjpvt, rcond, rank, 2469 nrr, ps, rcond, rank,
2467 work.fortran_vec (), lwork, prwork, 2470 work.fortran_vec (), lwork, prwork,
2468 info)); 2471 info));
2469 2472
2470 if (f77_exception_encountered) 2473 if (f77_exception_encountered)
2471 (*current_liboctave_error_handler) ("unrecoverable error in zgelsy"); 2474 (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
2472 else 2475 else
2473 { 2476 {
2474 lwork = static_cast<octave_idx_type> (std::real (work(0))); 2477 lwork = static_cast<octave_idx_type> (std::real (work(0)));
2475 work.resize (lwork); 2478 work.resize (lwork);
2476 2479
2477 F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, 2480 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
2478 nrr, pjpvt, rcond, rank, 2481 nrr, ps, rcond, rank,
2479 work.fortran_vec (), lwork, 2482 work.fortran_vec (), lwork,
2480 prwork, info)); 2483 prwork, info));
2481 2484
2482 if (f77_exception_encountered) 2485 if (f77_exception_encountered)
2483 (*current_liboctave_error_handler) 2486 (*current_liboctave_error_handler)
2484 ("unrecoverable error in zgelsy"); 2487 ("unrecoverable error in zgelss");
2485 else 2488 else
2486 { 2489 {
2487 retval.resize (n, nrhs); 2490 retval.resize (n, nrhs);
2488 for (octave_idx_type j = 0; j < nrhs; j++) 2491 for (octave_idx_type j = 0; j < nrhs; j++)
2489 for (octave_idx_type i = 0; i < n; i++) 2492 for (octave_idx_type i = 0; i < n; i++)
2558 for (octave_idx_type i = 0; i < m; i++) 2561 for (octave_idx_type i = 0; i < m; i++)
2559 result.elem (i) = b.elem (i); 2562 result.elem (i) = b.elem (i);
2560 2563
2561 Complex *presult = result.fortran_vec (); 2564 Complex *presult = result.fortran_vec ();
2562 2565
2563 Array<octave_idx_type> jpvt (n); 2566 octave_idx_type len_s = m < n ? m : n;
2564 octave_idx_type *pjpvt = jpvt.fortran_vec (); 2567 Array<double> s (len_s);
2568 double *ps = s.fortran_vec ();
2565 2569
2566 double rcond = -1.0; 2570 double rcond = -1.0;
2567 2571
2568 Array<double> rwork (2 * n); 2572 octave_idx_type lrwork = (5 * (m < n ? m : n)) - 4;
2573 lrwork = lrwork > 1 ? lrwork : 1;
2574 Array<double> rwork (lrwork);
2569 double *prwork = rwork.fortran_vec (); 2575 double *prwork = rwork.fortran_vec ();
2570 2576
2571 // Ask ZGELSY what the dimension of WORK should be. 2577 // Ask ZGELSS what the dimension of WORK should be.
2572 2578
2573 octave_idx_type lwork = -1; 2579 octave_idx_type lwork = -1;
2574 2580
2575 Array<Complex> work (1); 2581 Array<Complex> work (1);
2576 2582
2577 F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, 2583 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
2578 nrr, pjpvt, rcond, rank, 2584 nrr, ps, rcond, rank,
2579 work.fortran_vec (), lwork, prwork, 2585 work.fortran_vec (), lwork, prwork,
2580 info)); 2586 info));
2581 2587
2582 if (f77_exception_encountered) 2588 if (f77_exception_encountered)
2583 (*current_liboctave_error_handler) ("unrecoverable error in zgelsy"); 2589 (*current_liboctave_error_handler) ("unrecoverable error in zgelss");
2584 else 2590 else
2585 { 2591 {
2586 lwork = static_cast<int> (std::real (work(0))); 2592 lwork = static_cast<int> (std::real (work(0)));
2587 work.resize (lwork); 2593 work.resize (lwork);
2588 2594
2589 F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, 2595 F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult,
2590 nrr, pjpvt, rcond, rank, 2596 nrr, ps, rcond, rank,
2591 work.fortran_vec (), lwork, 2597 work.fortran_vec (), lwork,
2592 prwork, info)); 2598 prwork, info));
2593 2599
2594 if (f77_exception_encountered) 2600 if (f77_exception_encountered)
2595 (*current_liboctave_error_handler) 2601 (*current_liboctave_error_handler)
2596 ("unrecoverable error in zgelsy"); 2602 ("unrecoverable error in zgelss");
2597 else 2603 else
2598 { 2604 {
2599 retval.resize (n); 2605 retval.resize (n);
2600 for (octave_idx_type i = 0; i < n; i++) 2606 for (octave_idx_type i = 0; i < n; i++)
2601 retval.elem (i) = result.elem (i); 2607 retval.elem (i) = result.elem (i);