Mercurial > hg > octave-terminal
changeset 636:fae2bd91c027
[project @ 1994-08-23 18:39:50 by jwe]
author | jwe |
---|---|
date | Tue, 23 Aug 1994 18:39:50 +0000 |
parents | 5338832d2cf6 |
children | ae87e84f1be3 |
files | src/arith-ops.cc src/balance.cc src/chol.cc src/colloc.cc src/dassl.cc src/det.cc src/eig.cc src/expm.cc src/fft.cc src/file-io.cc src/find.cc src/fsolve.cc src/givens.cc src/hess.cc src/ifft.cc src/input.cc src/inv.cc src/log.cc src/lsode.cc src/lu.cc src/minmax.cc src/npsol.cc src/octave.cc src/pt-const.h src/pt-exp-base.cc src/qpsol.cc src/qr.cc src/quad.cc src/qzval.cc src/rand.cc src/schur.cc src/sort.cc src/svd.cc src/syl.cc src/sysdep.cc src/tc-rep.cc src/tc-rep.h src/utils.cc src/variables.cc |
diffstat | 39 files changed, 1728 insertions(+), 1770 deletions(-) [+] |
line wrap: on
line diff
--- a/src/arith-ops.cc +++ b/src/arith-ops.cc @@ -26,7 +26,6 @@ #endif #include <ctype.h> -#include <setjmp.h> #include <math.h> #include <Complex.h>
--- a/src/balance.cc +++ b/src/balance.cc @@ -39,6 +39,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -74,8 +75,9 @@ char *bal_job; int my_nargin; // # args w/o optional string arg - // determine if balancing option is listed - // set my_nargin to the number of matrix inputs +// Determine if balancing option is listed. Set my_nargin to the +// number of matrix inputs. + if (args(nargin-1).is_string ()) { bal_job = args(nargin-1).string_value (); @@ -87,30 +89,15 @@ my_nargin = nargin-1; } - tree_constant arg = args(1).make_numeric (); - int a_nr = arg.rows (); - int a_nc = arg.columns (); + tree_constant arg_a = args(1); + + int a_nr = arg_a.rows (); + int a_nc = arg_a.columns (); // Check argument 1 dimensions. - if (a_nr == 0 || a_nc == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - warning ("balance: argument is empty matrix"); - - Matrix m; - retval.resize (2); - retval(0) = m; - retval(1) = m; - } - else - error ("balance: empty matrix is invalid as argument"); - - return retval; - } + if (empty_arg ("balance", a_nr, a_nc) < 0) + return retval; if (a_nr != a_nc) { @@ -122,38 +109,23 @@ Matrix aa; ComplexMatrix caa; - if (arg.is_complex_type ()) - { - if (arg.is_matrix_type ()) - caa=arg.complex_matrix_value (); - else - { - caa.resize (1, 1); - caa.elem (0, 0) = arg.complex_value (); - } - } + if (arg_a.is_complex_type ()) + caa = arg_a.complex_matrix_value (); else - { - if (arg.is_matrix_type ()) - aa = arg.matrix_value (); - else - { - double d = arg.double_value (); - aa.resize (1, 1); - aa.elem (0, 0) = d; - } - } + aa = arg_a.matrix_value (); + + if (error_state) + return retval; // Treat AEP/ GEP cases. switch (my_nargin) { case 1: - + // Algebraic eigenvalue problem. - retval.resize (nargout ? nargout : 1); - if (arg.is_complex_type ()) + if (arg_a.is_complex_type ()) { ComplexAEPBALANCE result (caa, bal_job); @@ -161,8 +133,8 @@ retval(0) = result.balanced_matrix (); else { + retval(1) = result.balanced_matrix (); retval(0) = result.balancing_matrix (); - retval(1) = result.balanced_matrix (); } } else @@ -173,27 +145,26 @@ retval(0) = result.balanced_matrix (); else { + retval(1) = result.balanced_matrix (); retval(0) = result.balancing_matrix (); - retval(1) = result.balanced_matrix (); } } break; + case 2: - + { // Generalized eigenvalue problem. - { - retval.resize (nargout ? nargout : 1); - // 1st we have to check argument 2 dimensions and type... - tree_constant brg = args(2).make_numeric (); - int b_nr = brg.rows (); - int b_nc = brg.columns (); + tree_constant arg_b = args(2); + + int b_nr = arg_b.rows (); + int b_nc = arg_b.columns (); // Check argument 2 dimensions -- must match arg 1. - if ((b_nr != b_nc) || (b_nr != a_nr)) + if (b_nr != b_nc || b_nr != a_nr) { gripe_nonconformant (); return retval; @@ -204,42 +175,29 @@ Matrix bb; ComplexMatrix cbb; - if (brg.is_complex_type ()) - { - if (brg.is_matrix_type ()) - cbb = brg.complex_matrix_value (); - else - { - cbb.resize (1, 1); - cbb.elem (0, 0) = brg.complex_value (); - } - } + if (arg_b.is_complex_type ()) + cbb = arg_b.complex_matrix_value (); else - { - if (brg.is_matrix_type ()) - bb = brg.matrix_value (); - else - { - double d = brg.double_value (); - bb.resize (1, 1); - bb.elem (0, 0) = d; - } - } - + bb = arg_b.matrix_value (); + + if (error_state) + return retval; + // Both matrices loaded, now let's check what kind of arithmetic: - if (arg.is_complex_type () || brg.is_complex_type ()) + if (arg_a.is_complex_type () || arg_b.is_complex_type ()) { - if (arg.is_real_type ()) + if (arg_a.is_real_type ()) caa = aa; - else if (brg.is_real_type ()) + + if (arg_b.is_real_type ()) cbb = bb; // Compute magnitudes of elements for balancing purposes. // Surely there's a function I can call someplace! for (int i = 0; i < a_nr; i++) - for (int j = 0; j < a_nr; j++) + for (int j = 0; j < a_nc; j++) { aa.elem (i, j) = abs (caa.elem (i, j)); bb.elem (i, j) = abs (cbb.elem (i, j)); @@ -248,7 +206,7 @@ GEPBALANCE result (aa, bb, bal_job); - if (arg.is_complex_type () || brg.is_complex_type ()) + if (arg_a.is_complex_type () || arg_b.is_complex_type ()) { caa = result.left_balancing_matrix () * caa * result.right_balancing_matrix (); @@ -263,16 +221,19 @@ warning ("balance: should use two output arguments"); retval(0) = caa; break; + case 2: + retval(1) = cbb; retval(0) = caa; - retval(1) = cbb; break; + case 4: - retval(0) = result.left_balancing_matrix (); - retval(1) = result.right_balancing_matrix (); + retval(3) = cbb; retval(2) = caa; - retval(3) = cbb; + retval(1) = result.right_balancing_matrix (); + retval(0) = result.left_balancing_matrix (); break; + default: error ("balance: invalid number of output arguments"); break; @@ -287,16 +248,19 @@ warning ("balance: should use two output arguments"); retval(0) = result.balanced_a_matrix (); break; + case 2: + retval(1) = result.balanced_b_matrix (); retval(0) = result.balanced_a_matrix (); - retval(1) = result.balanced_b_matrix (); break; + case 4: - retval(0) = result.left_balancing_matrix (); - retval(1) = result.right_balancing_matrix (); + retval(3) = result.balanced_b_matrix (); retval(2) = result.balanced_a_matrix (); - retval(3) = result.balanced_b_matrix (); + retval(1) = result.right_balancing_matrix (); + retval(0) = result.left_balancing_matrix (); break; + default: error ("balance: invalid number of output arguments"); break; @@ -304,10 +268,12 @@ } } break; + default: error ("balance requires one (AEP) or two (GEP) numeric arguments"); break; } + return retval; }
--- a/src/chol.cc +++ b/src/chol.cc @@ -32,6 +32,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -40,68 +41,51 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2 || nargout > 1) + if (args.length () != 2 || nargout > 1) { print_usage ("chol"); return retval; } - tree_constant tmp = args(1).make_numeric (); + tree_constant arg = args(1); - int nr = tmp.rows (); - int nc = tmp.columns (); + int nr = arg.rows (); + int nc = arg.columns (); - if (nr == 0 || nc == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - gripe_empty_arg ("chol", 0); + if (empty_arg ("chol", nr, nc) < 0) + return retval; - retval.resize (1, Matrix ()); - } - else - gripe_empty_arg ("chol", 1); - - return retval; - } - - if (tmp.is_real_matrix ()) + if (arg.is_real_type ()) { - Matrix m = tmp.matrix_value (); - int info; - CHOL fact (m, info); - if (info != 0) - error ("chol: matrix not positive definite"); - else - retval = fact.chol_matrix (); + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int info; + CHOL fact (m, info); + if (info != 0) + error ("chol: matrix not positive definite"); + else + retval = fact.chol_matrix (); + } } - else if (tmp.is_complex_matrix ()) + else if (arg.is_complex_type ()) { - ComplexMatrix m = tmp.complex_matrix_value (); - int info; - ComplexCHOL fact (m, info); - if (info != 0) - error ("chol: matrix not positive definite"); - else - retval = fact.chol_matrix (); - } - else if (tmp.is_real_scalar ()) - { - double d = tmp.double_value (); - retval = d; - } - else if (tmp.is_complex_scalar ()) - { - Complex c = tmp.complex_value (); - retval = c; + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int info; + ComplexCHOL fact (m, info); + if (info != 0) + error ("chol: matrix not positive definite"); + else + retval = fact.chol_matrix (); + } } else { - gripe_wrong_type_arg ("chol", tmp); + gripe_wrong_type_arg ("chol", arg); } return retval;
--- a/src/colloc.cc +++ b/src/colloc.cc @@ -52,7 +52,12 @@ return retval; } - int ncol = NINT (args(1).double_value ()); + double tmp = args(1).double_value (); + + if (error_state) + return retval; + + int ncol = NINT (tmp); if (ncol < 0) { error ("colloc: first argument must be non-negative"); @@ -74,6 +79,7 @@ } char *s = args(i).string_value (); + if (s && (((*s == 'R' || *s == 'r') && strlen (s) == 1) || strcmp (s, "right") == 0)) { @@ -111,12 +117,10 @@ Matrix B = wts.second (); ColumnVector q = wts.quad_weights (); - retval.resize (4); - + retval(3) = q; + retval(2) = B; + retval(1) = A; retval(0) = r; - retval(1) = A; - retval(2) = B; - retval(3) = q; return retval; }
--- a/src/dassl.cc +++ b/src/dassl.cc @@ -95,7 +95,7 @@ { retval = tmp(0).vector_value (); - if (retval.length () == 0) + if (error_state || retval.length () == 0) gripe_user_supplied_eval ("dassl"); } else @@ -131,13 +131,41 @@ return retval; ColumnVector state = args(2).vector_value (); + + if (error_state) + { + error ("dassl: expecting state vector as second argument"); + return retval; + } + ColumnVector deriv = args(3).vector_value (); + + if (error_state) + { + error ("dassl: expecting derivative vector as third argument"); + return retval; + } + ColumnVector out_times = args(4).vector_value (); + + if (error_state) + { + error ("dassl: expecting output time vector as fourth argument"); + return retval; + } + ColumnVector crit_times; int crit_times_set = 0; if (nargin > 5) { crit_times = args(5).vector_value (); + + if (error_state) + { + error ("dassl: expecting critical time vector as fifth argument"); + return retval; + } + crit_times_set = 1; } @@ -280,20 +308,25 @@ if (nargin == 1) { print_dassl_option_list (); + return retval; } else if (nargin == 3) { - if (args(1).is_string ()) + char *keyword = args(1).string_value (); + + if (! error_state) { - char *keyword = args(1).string_value (); double val = args(2).double_value (); - do_dassl_option (keyword, val); + + if (! error_state) + { + do_dassl_option (keyword, val); + return retval; + } } - else - print_usage ("dassl_options"); } - else - print_usage ("dassl_options"); + + print_usage ("dassl_options"); return retval; }
--- a/src/det.cc +++ b/src/det.cc @@ -32,6 +32,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -48,31 +49,39 @@ return retval; } - tree_constant tmp = args(1).make_numeric ();; + tree_constant arg = args(1); - int nr = tmp.rows (); - int nc = tmp.columns (); - if (nr == 0 || nc == 0) + int nr = arg.rows (); + int nc = arg.columns (); + + if (nr == 0 && nc == 0) { - int flag = user_pref.propagate_empty_matrices; - if (flag < 0) - gripe_empty_arg ("det", 0); - else if (flag == 0) - gripe_empty_arg ("det", 1); + retval = 1.0; + return retval; } - if (nr == 0 && nc == 0) - return 1.0; + if (empty_arg ("det", nr, nc) < 0) + return retval; - if (tmp.is_real_matrix ()) + if (nr != nc) { - Matrix m = tmp.matrix_value (); - if (m.rows () == m.columns ()) + gripe_square_matrix_required ("det"); + return retval; + } + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) { int info; double rcond = 0.0; + DET det = m.determinant (info, rcond); + double d = 0.0; + if (info == -1) warning ("det: matrix singular to machine precision, rcond = %g", rcond); @@ -81,18 +90,20 @@ retval = d; } - else - gripe_square_matrix_required ("det"); } - else if (tmp.is_complex_matrix ()) + else if (arg.is_complex_matrix ()) { - ComplexMatrix m = tmp.complex_matrix_value (); - if (m.rows () == m.columns ()) + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) { int info; double rcond = 0.0; + ComplexDET det = m.determinant (info, rcond); + Complex c = 0.0; + if (info == -1) warning ("det: matrix singular to machine precision, rcond = %g", rcond); @@ -101,22 +112,10 @@ retval = c; } - else - gripe_square_matrix_required ("det"); - } - else if (tmp.is_real_scalar ()) - { - double d = tmp.double_value (); - retval = d; - } - else if (tmp.is_complex_scalar ()) - { - Complex c = tmp.complex_value (); - retval = c; } else { - gripe_wrong_type_arg ("det", tmp); + gripe_wrong_type_arg ("det", arg); } return retval;
--- a/src/eig.cc +++ b/src/eig.cc @@ -31,6 +31,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -39,37 +40,21 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2 || nargout > 2) + if (args.length () != 2 || nargout > 2) { print_usage ("eig"); return retval; } - tree_constant arg = args(1).make_numeric (); + tree_constant arg = args(1); - int a_nr = arg.rows (); - int a_nc = arg.columns (); + int nr = arg.rows (); + int nc = arg.columns (); - if (a_nr == 0 || a_nc == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - gripe_empty_arg ("eig", 0); - Matrix m; - retval(1) = m; - retval(0) = m; - } - else - gripe_empty_arg ("eig", 1); + if (empty_arg ("eig", nr, nc) < 0) + return retval; - return retval; - } - - if (a_nr != a_nc) + if (nr != nc) { gripe_square_matrix_required ("eig"); return retval; @@ -78,27 +63,24 @@ Matrix tmp; ComplexMatrix ctmp; EIG result; - if (arg.is_real_scalar ()) - { - tmp.resize (1, 1); - tmp.elem (0, 0) = arg.double_value (); - result = EIG (tmp); - } - else if (arg.is_real_matrix ()) + + if (arg.is_real_type ()) { tmp = arg.matrix_value (); - result = EIG (tmp); + + if (error_state) + return retval; + else + result = EIG (tmp); } - else if (arg.is_complex_scalar ()) - { - ctmp.resize (1, 1); - ctmp.elem (0, 0) = arg.complex_value (); - result = EIG (ctmp); - } - else if (arg.is_complex_matrix ()) + else if (arg.is_complex_type ()) { ctmp = arg.complex_matrix_value (); - result = EIG (ctmp); + + if (error_state) + return retval; + else + result = EIG (ctmp); } else {
--- a/src/expm.cc +++ b/src/expm.cc @@ -40,6 +40,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -65,7 +66,7 @@ return retval; } - tree_constant tmp = args(1).make_numeric (); + tree_constant arg = args(1); // Constants for matrix exponential calculation. @@ -81,231 +82,223 @@ 1.9270852604185938e-9, }; - if (tmp.is_empty ()) + int nr = arg.rows (); + int nc = arg.columns (); + + if (empty_arg ("expm", nr, nc) < 0) + return retval; + + if (nr != nc) { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - gripe_empty_arg ("expm", 0); - - retval.resize (1, Matrix ()); - } - else gripe_empty_arg ("expm", 1); + gripe_square_matrix_required ("expm"); + return retval; } - else if (tmp.rows () != tmp.columns ()) - gripe_square_matrix_required ("expm"); - else - { - int i, j; - int n_cols = tmp.columns (); + + int i, j; - char* balance_job = "B"; // variables for balancing + char* balance_job = "B"; // variables for balancing + + int sqpow; // power for scaling and squaring + double inf_norm; // norm of preconditioned matrix + int minus_one_j; // used in computing pade approx - int sqpow; // power for scaling and squaring - double inf_norm; // norm of preconditioned matrix - int minus_one_j; // used in computing pade approx + if (arg.is_real_type ()) + { + +// Compute the exponential. - if (tmp.is_complex_matrix ()) - { - ComplexMatrix m = tmp.complex_matrix_value (); - Complex trshift = 0.0; // trace shift value + Matrix m = arg.matrix_value (); + + if (error_state) + return retval; + + double trshift = 0; // trace shift value // Preconditioning step 1: trace normalization. - for (i = 0; i < n_cols; i++) - trshift += m.elem (i, i); - trshift /= n_cols; - for (i = 0; i < n_cols; i++) - m.elem (i, i) -= trshift; + for (i = 0; i < nc; i++) + trshift += m.elem (i, i); + trshift /= nc; + for (i = 0; i < nc; i++) + m.elem (i, i) -= trshift; -// Preconditioning step 2: eigenvalue balancing. +// Preconditioning step 2: balancing. - ComplexAEPBALANCE mbal (m, balance_job); - m = mbal.balanced_matrix (); - ComplexMatrix d = mbal.balancing_matrix (); + AEPBALANCE mbal (m, balance_job); + m = mbal.balanced_matrix (); + Matrix d = mbal.balancing_matrix (); // Preconditioning step 3: scaling. - ColumnVector work (n_cols); - inf_norm = F77_FCN (zlange) ("I", &n_cols, &n_cols, m. - fortran_vec (), &n_cols, - work.fortran_vec ()); + ColumnVector work(nc); + inf_norm = F77_FCN (dlange) ("I", &nc, &nc, + m.fortran_vec (), &nc, + work.fortran_vec ()); - sqpow = (int) (1.0 + log (inf_norm) / log (2.0)); + sqpow = (int) (1.0 + log (inf_norm) / log (2.0)); // Check whether we need to square at all. - if (sqpow < 0) - sqpow = 0; - else - { - for (inf_norm = 1.0, i = 0; i < sqpow; i++) - inf_norm *= 2.0; + if (sqpow < 0) + sqpow = 0; + else + { + for (inf_norm = 1.0, i = 0; i < sqpow; i++) + inf_norm *= 2.0; - m = m / inf_norm; - } + m = m / inf_norm; + } // npp, dpp: pade' approx polynomial matrices. - ComplexMatrix npp (n_cols, n_cols, 0.0); - ComplexMatrix dpp = npp; + Matrix npp (nc, nc, 0.0); + Matrix dpp = npp; + +// now powers a^8 ... a^1. + + minus_one_j = -1; + for (j = 7; j >= 0; j--) + { + npp = m * npp + m * padec[j]; + dpp = m * dpp + m * (minus_one_j * padec[j]); + minus_one_j *= -1; + } +// Zero power. + + dpp = -dpp; + for(j = 0; j < nc; j++) + { + npp.elem (j, j) += 1.0; + dpp.elem (j, j) += 1.0; + } + +// Compute pade approximation = inverse (dpp) * npp. + + Matrix result = dpp.solve (npp); + +// Reverse preconditioning step 3: repeated squaring. + + while (sqpow) + { + result = result * result; + sqpow--; + } + +// Reverse preconditioning step 2: inverse balancing. + + result = result.transpose(); + d = d.transpose (); + result = result * d; + result = d.solve (result); + result = result.transpose (); + +// Reverse preconditioning step 1: fix trace normalization. + + result = result * exp (trshift); + + retval = result; + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (error_state) + return retval; + + Complex trshift = 0.0; // trace shift value + +// Preconditioning step 1: trace normalization. + + for (i = 0; i < nc; i++) + trshift += m.elem (i, i); + trshift /= nc; + for (i = 0; i < nc; i++) + m.elem (i, i) -= trshift; + +// Preconditioning step 2: eigenvalue balancing. + + ComplexAEPBALANCE mbal (m, balance_job); + m = mbal.balanced_matrix (); + ComplexMatrix d = mbal.balancing_matrix (); + +// Preconditioning step 3: scaling. + + ColumnVector work (nc); + inf_norm = F77_FCN (zlange) ("I", &nc, &nc, m. + fortran_vec (), &nc, + work.fortran_vec ()); + + sqpow = (int) (1.0 + log (inf_norm) / log (2.0)); + +// Check whether we need to square at all. + + if (sqpow < 0) + sqpow = 0; + else + { + for (inf_norm = 1.0, i = 0; i < sqpow; i++) + inf_norm *= 2.0; + + m = m / inf_norm; + } + +// npp, dpp: pade' approx polynomial matrices. + + ComplexMatrix npp (nc, nc, 0.0); + ComplexMatrix dpp = npp; // Now powers a^8 ... a^1. - minus_one_j = -1; - for (j = 7; j >= 0; j--) - { - npp = m * npp + m * padec[j]; - dpp = m * dpp + m * (minus_one_j * padec[j]); - minus_one_j *= -1; - } + minus_one_j = -1; + for (j = 7; j >= 0; j--) + { + npp = m * npp + m * padec[j]; + dpp = m * dpp + m * (minus_one_j * padec[j]); + minus_one_j *= -1; + } // Zero power. - dpp = -dpp; - for (j = 0; j < n_cols; j++) - { - npp.elem (j, j) += 1.0; - dpp.elem (j, j) += 1.0; - } + dpp = -dpp; + for (j = 0; j < nc; j++) + { + npp.elem (j, j) += 1.0; + dpp.elem (j, j) += 1.0; + } // Compute pade approximation = inverse (dpp) * npp. - ComplexMatrix result = dpp.solve (npp); + ComplexMatrix result = dpp.solve (npp); // Reverse preconditioning step 3: repeated squaring. - while (sqpow) - { - result = result * result; - sqpow--; - } + while (sqpow) + { + result = result * result; + sqpow--; + } // reverse preconditioning step 2: inverse balancing XXX FIXME XXX: // should probably do this with lapack calls instead of a complete // matrix inversion. - result = result.transpose (); - d = d.transpose (); - result = result * d; - result = d.solve (result); - result = result.transpose (); + result = result.transpose (); + d = d.transpose (); + result = result * d; + result = d.solve (result); + result = result.transpose (); // Reverse preconditioning step 1: fix trace normalization. - result = result * exp (trshift); - - retval = result; - } - else if (tmp.is_complex_scalar ()) - { - Complex c = tmp.complex_value (); - retval = exp (c); - } - else if (tmp.is_real_matrix ()) - { - -// Compute the exponential. - - Matrix m = tmp.matrix_value (); - - double trshift = 0; // trace shift value - -// Preconditioning step 1: trace normalization. - - for (i = 0; i < n_cols; i++) - trshift += m.elem (i, i); - trshift /= n_cols; - for (i = 0; i < n_cols; i++) - m.elem (i, i) -= trshift; - -// Preconditioning step 2: balancing. - - AEPBALANCE mbal (m, balance_job); - m = mbal.balanced_matrix (); - Matrix d = mbal.balancing_matrix (); - -// Preconditioning step 3: scaling. - - ColumnVector work(n_cols); - inf_norm = F77_FCN (dlange) ("I", &n_cols, &n_cols, - m.fortran_vec (), &n_cols, - work.fortran_vec ()); - - sqpow = (int) (1.0 + log (inf_norm) / log (2.0)); - -// Check whether we need to square at all. - - if (sqpow < 0) - sqpow = 0; - else - { - for (inf_norm = 1.0, i = 0; i < sqpow; i++) - inf_norm *= 2.0; - - m = m / inf_norm; - } - -// npp, dpp: pade' approx polynomial matrices. + result = result * exp (trshift); - Matrix npp (n_cols, n_cols, 0.0); - Matrix dpp = npp; - -// now powers a^8 ... a^1. - - minus_one_j = -1; - for (j = 7; j >= 0; j--) - { - npp = m * npp + m * padec[j]; - dpp = m * dpp + m * (minus_one_j * padec[j]); - minus_one_j *= -1; - } -// Zero power. - - dpp = -dpp; - for(j = 0; j < n_cols; j++) - { - npp.elem (j, j) += 1.0; - dpp.elem (j, j) += 1.0; - } - -// Compute pade approximation = inverse (dpp) * npp. - - Matrix result = dpp.solve (npp); - -// Reverse preconditioning step 3: repeated squaring. + retval = result; + } + else + { + gripe_wrong_type_arg ("expm", arg); + } - while (sqpow) - { - result = result * result; - sqpow--; - } - -// Reverse preconditioning step 2: inverse balancing. - - result = result.transpose(); - d = d.transpose (); - result = result * d; - result = d.solve (result); - result = result.transpose (); - -// Reverse preconditioning step 1: fix trace normalization. - - result = result * exp (trshift); - - retval = result; - } - else if (tmp.is_real_scalar ()) - { - double d = tmp.double_value (); - retval = exp (d); - } - else - { - gripe_wrong_type_arg ("expm", tmp); - } - } return retval; }
--- a/src/fft.cc +++ b/src/fft.cc @@ -32,6 +32,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -40,51 +41,40 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2) + if (args.length () != 2) { print_usage ("fft"); return retval; } - tree_constant tmp = args(1).make_numeric ();; - - if (tmp.rows () == 0 || tmp.columns () == 0) + tree_constant arg = args(1); + + if (empty_arg ("fft", arg.rows (), arg.columns ()) < 0) + return retval; + + if (arg.is_real_type ()) { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) + Matrix m = arg.matrix_value (); + + if (! error_state) { - if (flag < 0) - gripe_empty_arg ("fft", 0); - - retval.resize (1, Matrix ()); + ComplexMatrix mfft = m.fourier (); + retval = mfft; } - else - gripe_empty_arg ("fft", 1); - - return retval; } - - if (tmp.is_real_matrix ()) + else if (arg.is_complex_type ()) { - Matrix m = tmp.matrix_value (); - ComplexMatrix mfft = m.fourier (); - retval = mfft; - } - else if (tmp.is_complex_matrix ()) - { - ComplexMatrix m = tmp.complex_matrix_value (); - ComplexMatrix mfft = m.fourier (); - retval = mfft; - } - else if (tmp.is_scalar_type ()) - { - error ("fft: invalid scalar argument"); + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix mfft = m.fourier (); + retval = mfft; + } } else { - gripe_wrong_type_arg ("fft", tmp); + gripe_wrong_type_arg ("fft", arg); } return retval;
--- a/src/file-io.cc +++ b/src/file-io.cc @@ -190,41 +190,43 @@ file_list.next (p); } } - else if (arg.is_scalar_type ()) + else { double file_num = arg.double_value (); - if ((double) NINT (file_num) != file_num) - error ("file number not an integer value"); - else + + if (! error_state) { - Pix p = file_list.first (); - file_info file; - for (int i = 0; i < file_count; i++) + if ((double) NINT (file_num) != file_num) + error ("file number not an integer value"); + else { - file = file_list (p); - if (file.number () == file_num) - return p; - file_list.next (p); + Pix p = file_list.first (); + file_info file; + for (int i = 0; i < file_count; i++) + { + file = file_list (p); + if (file.number () == file_num) + return p; + file_list.next (p); + } + error ("no file with that number"); } - error ("no file with that number"); } - } - else - error ("inapproriate file specifier"); + else + error ("inapproriate file specifier"); + } return 0; } static Pix -fopen_file_for_user (const tree_constant& arg, const char *mode, +fopen_file_for_user (const char *name, const char *mode, const char *warn_for) { - char *file_name = arg.string_value (); - - FILE *file_ptr = fopen (file_name, mode); + FILE *file_ptr = fopen (name, mode); if (file_ptr) { - file_info file (++file_count, file_name, file_ptr, mode); + file_info file (++file_count, name, file_ptr, mode); file_list.append (file); Pix p = file_list.first (); @@ -233,19 +235,19 @@ for (int i = 0; i < file_count; i++) { file_from_list = file_list (p); - if (strcmp (file_from_list.name (), file_name) == 0) + if (strcmp (file_from_list.name (), name) == 0) return p; file_list.next (p); } } - error ("%s: unable to open file `%s'", warn_for, file_name); + error ("%s: unable to open file `%s'", warn_for, name); return 0; } static Pix -file_io_get_file (const tree_constant arg, const char *mode, +file_io_get_file (const tree_constant& arg, const char *mode, const char *warn_for) { Pix p = return_valid_file (arg); @@ -262,12 +264,12 @@ if (status == 0) { if ((buffer.st_mode & S_IFREG) == S_IFREG) - p = fopen_file_for_user (arg, mode, warn_for); + p = fopen_file_for_user (name, mode, warn_for); else error ("%s: invalid file type", warn_for); } else if (status < 0 && *mode != 'r') - p = fopen_file_for_user (arg, mode, warn_for); + p = fopen_file_for_user (name, mode, warn_for); else error ("%s: can't stat file `%s'", warn_for, name); } @@ -315,13 +317,12 @@ file_list.del (p); file_count--; - retval.resize (1); if (success == 0) - retval(0) = tree_constant (1.0); // succeeded + retval(0) = 1.0; // succeeded else { error ("fclose: error on closing file"); - retval(0) = tree_constant (0.0); // failed + retval(0) = 0.0; // failed } return retval; @@ -367,13 +368,12 @@ else success = fflush (file.fptr ()); - retval.resize (1); if (success == 0) - retval(0) = tree_constant (1.0); // succeeded + retval(0) = 1.0; // succeeded else { error ("fflush: write error"); - retval(0) = tree_constant (0.0); // failed + retval(0) = 0.0; // failed } return retval; @@ -421,38 +421,35 @@ if (! p) return retval; - int length = 0; - if (args(2).is_scalar_type ()) + + double dlen = args(2).double_value (); + + if (error_state) + return retval; + + int length = NINT (dlen); + + if ((double) length != dlen) { - length = (int) args(2).double_value (); - if ((double) NINT (length) != length) - { - error ("fgets: length not an integer value"); - return retval; - } + error ("fgets: length not an integer value"); + return retval; } file_info file = file_list (p); - char string[length+1]; + char string [length+1]; char *success = fgets (string, length+1, file.fptr ()); if (! success) { - retval.resize (1); - retval(0) = tree_constant (-1.0); + retval(0) = -1.0; return retval; } if (nargout == 2) - { - retval.resize (2); - retval(1) = tree_constant ((double) strlen (string)); - } - else - retval.resize (1); + retval(1) = (double) strlen (string); - retval(0) = tree_constant (string); + retval(0) = string; return retval; } @@ -500,8 +497,7 @@ { file_info file = file_list (p); - retval.resize (1); - retval(0) = tree_constant ((double) file.number ()); + retval(0) = (double) file.number (); return retval; } @@ -541,8 +537,7 @@ file_info file (number, name, file_ptr, mode); file_list.append (file); - retval.resize (1); - retval(0) = tree_constant ((double) number); + retval(0) = (double) number; return retval; } @@ -646,44 +641,55 @@ return retval; long origin = SEEK_SET; - long offset = 0; - if (args(2).is_scalar_type ()) + + double doff = args(2).double_value (); + + if (error_state) + return retval; + + long offset = NINT (doff); + + if ((double) offset != doff) { - offset = (long) args(2).double_value (); - if ((double) NINT (offset) != offset) + error ("fseek: offset not an integer value"); + return retval; + } + + if (nargin == 4) + { + double dorig = args(3).double_value (); + + if (error_state) + return retval; + + origin = NINT (dorig); + + if ((double) dorig != origin) { - error ("fseek: offset not an integer value"); + error ("fseek: origin not an integer value"); return retval; } - } - if (nargin == 4 && args(3).is_scalar_type ()) - { - origin = (long) args(3).double_value (); if (origin == -1) origin = SEEK_CUR; else if (origin == -2) origin = SEEK_END; else { - if ((double) NINT (origin) != origin) - { - error ("fseek: origin not an integer value"); - return retval; - } + error ("fseek: invalid value for origin"); + return retval; } } file_info file = file_list (p); int success = fseek (file.fptr (), offset, origin); - retval.resize (1); if (success == 0) - retval(0) = tree_constant (1.0); // succeeded + retval(0) = 1.0; // succeeded else { error ("fseek: file error"); - retval(0) = tree_constant (0.0); // failed + retval(0) = 0.0; // failed } return retval; @@ -718,8 +724,8 @@ { file_info file = file_list (p); long offset = ftell (file.fptr ()); - retval.resize (1); - retval(0) = tree_constant ((double) offset); + + retval(0) = (double) offset; if (offset == -1L) error ("ftell: write error"); @@ -781,13 +787,15 @@ return -1; } - if (! args(fmt_arg_count).is_scalar_type ()) + double tmp_len = args(fmt_arg_count++).double_value (); + + if (error_state) { error ("%s: `*' must be replaced by an integer", type); return -1; } - fmt << NINT (args(fmt_arg_count++).double_value ()); + fmt << NINT (tmp_len); s++; chars_from_fmt_str++; } @@ -820,13 +828,15 @@ return -1; } - if (! args(fmt_arg_count).is_scalar_type ()) + double tmp_len = args(fmt_arg_count++).double_value (); + + if (error_state) { error ("%s: `*' must be replaced by an integer", type); return -1; } - fmt << NINT (args(fmt_arg_count++).double_value ()); + fmt << NINT (tmp_len); s++; chars_from_fmt_str++; } @@ -860,72 +870,74 @@ switch (*s) { case 'd': case 'i': case 'o': case 'u': case 'x': case 'X': + { + double d = args(fmt_arg_count++).double_value (); - if (! args(fmt_arg_count).is_scalar_type ()) - goto invalid_conversion; - else - { - chars_from_fmt_str++; - fmt << *s << ends; - double d = args(fmt_arg_count++).double_value (); - if ((int) d != d) - goto invalid_conversion; - else - { - char *s = fmt.str (); - sb.form (s, (int) d); - delete [] s; - return chars_from_fmt_str; - } - } + int val = NINT (d); + + if (error_state || (double) val != d) + goto invalid_conversion; + else + { + chars_from_fmt_str++; + fmt << *s << ends; + char *tmp_fmt = fmt.str (); + sb.form (tmp_fmt, val); + delete [] tmp_fmt; + return chars_from_fmt_str; + } + } case 'e': case 'E': case 'f': case 'g': case 'G': + { + double val = args(fmt_arg_count++).double_value (); - if (! args(fmt_arg_count).is_scalar_type ()) - goto invalid_conversion; - else - { - chars_from_fmt_str++; - fmt << *s << ends; - char *s = fmt.str (); - sb.form (s, args(fmt_arg_count++).double_value ()); - delete [] s; - return chars_from_fmt_str; - } + if (error_state) + goto invalid_conversion; + else + { + chars_from_fmt_str++; + fmt << *s << ends; + char *tmp_fmt = fmt.str (); + sb.form (tmp_fmt, val); + delete [] tmp_fmt; + return chars_from_fmt_str; + } + } case 's': + { + char *val = args(fmt_arg_count++).string_value (); - if (! args(fmt_arg_count).is_string ()) - goto invalid_conversion; - else - { - chars_from_fmt_str++; - fmt << *s << ends; - char *s = fmt.str (); - sb.form (s, args(fmt_arg_count++).string_value ()); - delete [] s; - return chars_from_fmt_str; - } + if (error_state) + goto invalid_conversion; + else + { + chars_from_fmt_str++; + fmt << *s << ends; + char *tmp_fmt = fmt.str (); + sb.form (tmp_fmt, val); + delete [] tmp_fmt; + return chars_from_fmt_str; + } + } case 'c': + { + char *val = args(fmt_arg_count++).string_value (); - if (! args(fmt_arg_count).is_string ()) - goto invalid_conversion; - else - { - chars_from_fmt_str++; - fmt << *s << ends; - char *str = args(fmt_arg_count++).string_value (); - if (strlen (str) != 1) - goto invalid_conversion; - else - { - char *s = fmt.str (); - sb.form (s, *str); - delete [] s; - return chars_from_fmt_str; - } - } + if (error_state || strlen (val) != 1) + goto invalid_conversion; + else + { + chars_from_fmt_str++; + fmt << *s << ends; + char *tmp_fmt = fmt.str (); + sb.form (tmp_fmt, *val); + delete [] tmp_fmt; + return chars_from_fmt_str; + } + } default: goto invalid_format; @@ -1005,17 +1017,6 @@ if (strcmp (type, "fprintf") == 0) { - if (args(2).is_string ()) - { - fmt = args(2).string_value (); - fmt_arg_count++; - } - else - { - error ("%s: format must be a string", type); - return retval; - } - Pix p = file_io_get_file (args(1), "a+", type); if (! p) @@ -1031,17 +1032,25 @@ fmt = args(2).string_value (); - fmt_arg_count++; - } - else if (args(1).is_string ()) - { - fmt = args(1).string_value (); - fmt_arg_count++; + if (error_state) + { + error ("%s: format must be a string", type); + return retval; + } + + fmt_arg_count += 2; } else { - error ("%s: invalid format string", type); - return retval; + fmt = args(1).string_value (); + + if (error_state) + { + error ("%s: invalid format string", type); + return retval; + } + + fmt_arg_count++; } // Scan fmt for % escapes and print out the arguments. @@ -1093,9 +1102,8 @@ } else if (strcmp (type, "sprintf") == 0) { - retval.resize (1); char *msg = output_buf.str (); - retval(0) = tree_constant (msg); + retval(0) = msg; delete [] msg; } @@ -1165,7 +1173,7 @@ success = fscanf (fptr, str, &temp); delete [] str; if (success > 0 && store_value) - values(fmt_arg_count++) = tree_constant ((double) temp); + values(fmt_arg_count++) = (double) temp; } break; @@ -1178,7 +1186,7 @@ success = fscanf (fptr, str, &temp); delete [] str; if (success > 0 && store_value) - values(fmt_arg_count++) = tree_constant (temp); + values(fmt_arg_count++) = temp; } break; @@ -1197,8 +1205,7 @@ int c; - while ((c = getc (fptr)) != EOF - && (c == ' ' || c == '\n' || c != '\t')) + while ((c = getc (fptr)) != EOF && isspace (c)) ; // Don't count leading whitespace. if (c != EOF) @@ -1207,7 +1214,7 @@ for (;;) { c = getc (fptr); - if (c != EOF && c != ' ' && c != '\n' && c != '\t') + if (c != EOF && ! isspace (c)) string_width++; else break; @@ -1216,13 +1223,14 @@ fseek (fptr, original_position, SEEK_SET); } chars_from_fmt_str++; - char temp[string_width+1]; + char temp [string_width+1]; fmt << *s << ends; char *str = fmt.str (); success = fscanf (fptr, str, temp); delete [] str; + temp[string_width] = '\0'; if (success > 0 && store_value) - values(fmt_arg_count++) = tree_constant (temp); + values(fmt_arg_count++) = temp; } break; @@ -1231,7 +1239,7 @@ if (string_width < 1) string_width = 1; chars_from_fmt_str++; - char temp[string_width+1]; + char temp [string_width+1]; memset (temp, '\0', string_width+1); fmt << *s << ends; char *str = fmt.str (); @@ -1239,7 +1247,7 @@ delete [] str; temp[string_width] = '\0'; if (success > 0 && store_value) - values(fmt_arg_count++) = tree_constant (temp); + values(fmt_arg_count++) = temp; } break; @@ -1337,9 +1345,9 @@ if (strcmp (type, "scanf") != 0) { - if (args(2).is_string ()) - scanf_fmt = args(2).string_value (); - else + scanf_fmt = args(2).string_value (); + + if (error_state) { error ("%s: format must be a string", type); return retval; @@ -1417,7 +1425,7 @@ // Scan scanf_fmt for % escapes and assign the arguments. - retval.resize (nargout ? nargout : 1); + retval.resize (nargout); char *ptr = scanf_fmt; @@ -1568,9 +1576,9 @@ char *prec = "uchar"; if (nargin > 3) { - if (args(3).is_string ()) - prec = args(3).string_value (); - else + prec = args(3).string_value (); + + if (error_state) { error ("fread: precision must be a specified as a string"); return retval; @@ -1595,24 +1603,25 @@ { if (args(2).is_scalar_type ()) { - tree_constant tmpa = args(2).make_numeric (); - dnr = tmpa.double_value (); + dnr = args(2).double_value (); + + if (error_state) + return retval; + dnc = 1.0; } - else if (args(2).is_matrix_type ()) + else { ColumnVector tmp = args(2).vector_value (); - if (tmp.length () == 2) - { - dnr = tmp.elem (0); - dnc = tmp.elem (1); - } - else + if (error_state || tmp.length () != 2) { error ("fread: invalid size specification\n"); return retval; } + + dnr = tmp.elem (0); + dnc = tmp.elem (1); } if ((xisinf (dnr)) && (xisinf (dnc))) @@ -1660,14 +1669,9 @@ int count = m.read (fptr, prec); if (nargout > 1) - { - retval.resize (2); - retval(1) = tree_constant ((double) count); - } - else - retval.resize (1); + retval(1) = (double) count; - retval(0) = tree_constant (m); + retval(0) = m; return retval; } @@ -1730,9 +1734,9 @@ char *prec = "uchar"; if (nargin > 3) { - if (args(3).is_string ()) - prec = args(3).string_value (); - else + prec = args(3).string_value (); + + if (error_state) { error ("fwrite: precision must be a specified as a string"); return retval; @@ -1743,10 +1747,12 @@ Matrix m = args(2).matrix_value (); - int count = m.write (file.fptr (), prec); + if (! error_state) + { + int count = m.write (file.fptr (), prec); - retval.resize (1); - retval(0) = tree_constant ((double) count); + retval(0) = (double) count; + } return retval; } @@ -1790,7 +1796,6 @@ file_info file = file_list (p); - retval.resize (1); retval(0) = (double) feof (file.fptr ()); return retval; @@ -1839,14 +1844,9 @@ int ierr = ferror (file.fptr ()); if (nargout > 1) - { - retval.resize (2); - retval(1) = tree_constant ((double) ierr); - } - else - retval.resize (1); + retval(1) = (double) ierr; - retval(0) = tree_constant (strsave (strerror (ierr))); + retval(0) = strsave (strerror (ierr)); return retval; }
--- a/src/find.cc +++ b/src/find.cc @@ -35,10 +35,10 @@ const tree_constant& val, int nr, int nc, int nargout) { Octave_object retval; - retval.resize (nargout); switch (nargout) { + case 0: case 1: { int count = i_idx.length (); @@ -51,15 +51,19 @@ // retval(0) = tree_constant (tmp, (nr != 1)); } break; + case 3: retval(2) = val; +// Fall through! + case 2: + retval(1) = tree_constant (j_idx, 1); retval(0) = tree_constant (i_idx, 1); // If you want this to work more like Matlab, use the following line // instead of the previous one. // retval(0) = tree_constant (i_idx, (nr != 1)); - retval(1) = tree_constant (j_idx, 1); break; + default: panic_impossible (); break; @@ -81,8 +85,7 @@ if (m.elem (i, j) != 0.0) count++; - Matrix result; - Octave_object retval (nargout, result); + Octave_object retval (((nargout == 0) ? 1 : nargout), Matrix ()); if (count == 0) return retval; @@ -122,8 +125,7 @@ if (m.elem (i, j) != 0.0) count++; - Matrix result; - Octave_object retval (nargout, result); + Octave_object retval (((nargout == 0) ? 1 : nargout), Matrix ()); if (count == 0) return retval; @@ -163,49 +165,25 @@ return retval; } - nargout = (nargout == 0) ? 1 : nargout; - - retval.resize (nargout, Matrix ()); - - tree_constant tmp = args(1).make_numeric (); + tree_constant arg = args(1); - if (tmp.is_real_matrix ()) - { - Matrix m = tmp.matrix_value (); - return find_nonzero_elem_idx (m, nargout); - } - else if (tmp.is_real_scalar ()) + if (arg.is_real_type ()) { - double d = tmp.double_value (); - if (d != 0.0) - { - retval(0) = 1.0; - if (nargout > 1) - retval(1) = 1.0; - if (nargout > 2) - retval(2) = d; - } + Matrix m = arg.matrix_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (m, nargout); } - else if (tmp.is_complex_matrix ()) + else if (arg.is_complex_type ()) { - ComplexMatrix m = tmp.complex_matrix_value (); - return find_nonzero_elem_idx (m, nargout); - } - else if (tmp.is_complex_scalar ()) - { - Complex c = tmp.complex_value (); - if (c != 0.0) - { - retval(0) = 1.0; - if (nargout > 1) - retval(1) = 1.0; - if (nargout > 2) - retval(2) = c; - } + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (m, nargout); } else { - gripe_wrong_type_arg ("find", tmp); + gripe_wrong_type_arg ("find", arg); } return retval;
--- a/src/fsolve.cc +++ b/src/fsolve.cc @@ -51,23 +51,29 @@ case -1: info = -2; break; + case 0: info = -1; break; + case 1: break; + case 2: info = 4; break; + case 3: case 4: case 5: info = 3; break; + default: panic_impossible (); break; } + return info; } @@ -105,7 +111,7 @@ { retval = tmp(0).vector_value (); - if (retval.length () <= 0) + if (error_state || retval.length () <= 0) gripe_user_supplied_eval ("fsolve"); } else @@ -143,6 +149,12 @@ ColumnVector x = args(2).vector_value (); + if (error_state) + { + error ("fsolve: expecting vector as second argument"); + return retval; + } + if (nargin > 3) warning ("fsolve: ignoring extra arguments"); @@ -262,20 +274,25 @@ if (nargin == 1) { print_fsolve_option_list (); + return retval; } else if (nargin == 3) { - if (args(1).is_string ()) + char *keyword = args(1).string_value (); + + if (! error_state) { - char *keyword = args(1).string_value (); double val = args(2).double_value (); - do_fsolve_option (keyword, val); + + if (! error_state) + { + do_fsolve_option (keyword, val); + return retval; + } } - else - print_usage ("fsolve_options"); } - else - print_usage ("fsolve_options"); + + print_usage ("fsolve_options"); return retval; }
--- a/src/givens.cc +++ b/src/givens.cc @@ -47,16 +47,6 @@ Complex*, Complex*); } -// These aren't used? -#if 0 -int F77_FCN (dorgqr) (const int*, const int*, const int*, double*, - const int*, double*, double*, const int*, int*); - -int F77_FCN (zunghr) (const int*, const int*, const int*, Complex*, - const int*, Complex*, Complex*, const int*, - int*, long, long); -#endif - DEFUN_DLD ("givens", Fgivens, Sgivens, 3, 2, "G = givens (X, Y)\n\ \n\ @@ -67,108 +57,121 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 3 || nargout > 2) + if (args.length () != 3 || nargout > 2) { print_usage ("givens"); return retval; } - tree_constant arga = args(1).make_numeric (); - tree_constant argb = args(2).make_numeric (); + tree_constant arg_a = args(1); + tree_constant arg_b = args(2); + + if (! arg_a.is_scalar_type () && arg_b.is_scalar_type ()) + { + error("givens: requires two scalar arguments"); + return retval; + } + + Complex cx, cy; + double x, y; + + if (arg_a.is_complex_type ()) + { + cx = arg_a.complex_value (); - if (! arga.is_scalar_type () && argb.is_scalar_type ()) + if (error_state) + return retval; + } + else { - error("givens: requires two scalar arguments"); + x = arg_a.double_value (); + + if (error_state) + return retval; + + cx = x; // copy to complex just in case + } + + if (arg_b.is_complex_type ()) + { + cy = arg_b.complex_value (); + + if (error_state) + return retval; } else { - retval.resize (nargout ? nargout : 1); - - Complex cx, cy; - double x, y; + y = arg_b.double_value (); - if (arga.is_complex_type ()) - cx = arga.complex_value (); - else - { - x = arga.double_value (); - cx = x; // copy to complex just in case - } + if (error_state) + return retval; - if (argb.is_complex_type ()) - cy = argb.complex_value (); - else - { - y = argb.double_value (); - cy = y; // copy to complex just in case - } + cy = y; // copy to complex just in case + } // Now compute the rotation. - double cc; - if (arga.is_complex_type () || argb.is_complex_type ()) - { - Complex cs, temp_r; + double cc; + if (arg_a.is_complex_type () || arg_b.is_complex_type ()) + { + Complex cs, temp_r; - F77_FCN (zlartg) (&cx, &cy, &cc, &cs, &temp_r); + F77_FCN (zlartg) (&cx, &cy, &cc, &cs, &temp_r); - switch (nargout) - { - case 0: // output a matrix - case 1: - { - ComplexMatrix g (2, 2); - g.elem (0, 0) = cc; - g.elem (1, 1) = cc; - g.elem (0, 1) = cs; - g.elem (1, 0) = -conj (cs); + switch (nargout) + { + case 0: // output a matrix + case 1: + { + ComplexMatrix g (2, 2); + g.elem (0, 0) = cc; + g.elem (1, 1) = cc; + g.elem (0, 1) = cs; + g.elem (1, 0) = -conj (cs); - retval(0) = g; - } - break; + retval(0) = g; + } + break; - case 2: // output scalar values - retval(0) = tree_constant(cc); - retval(1) = tree_constant(cs); - break; + case 2: // output scalar values + retval(0) = cc; + retval(1) = cs; + break; - default: - error ("givens: invalid number of output arguments"); - break; - } + default: + error ("givens: invalid number of output arguments"); + break; } - else - { - double s, temp_r; - - F77_FCN (dlartg) (&x, &y, &cc, &s, &temp_r); + } + else + { + double s, temp_r; + + F77_FCN (dlartg) (&x, &y, &cc, &s, &temp_r); - switch (nargout) - { - case 0: // output a matrix - case 1: - { - Matrix g (2, 2); - g.elem (0, 0) = cc; - g.elem (1, 1) = cc; - g.elem (0, 1) = s; - g.elem (1, 0) = -s; + switch (nargout) + { + case 0: // output a matrix + case 1: + { + Matrix g (2, 2); + g.elem (0, 0) = cc; + g.elem (1, 1) = cc; + g.elem (0, 1) = s; + g.elem (1, 0) = -s; - retval(0) = g; - } - break; + retval(0) = g; + } + break; - case 2: // output scalar values - retval(0) = cc; - retval(1) = s; - break; + case 2: // output scalar values + retval(0) = cc; + retval(1) = s; + break; - default: - error ("givens: invalid number of output arguments"); - break; - } + default: + error ("givens: invalid number of output arguments"); + break; } }
--- a/src/hess.cc +++ b/src/hess.cc @@ -32,6 +32,7 @@ #include "user-prefs.h" #include "error.h" #include "gripes.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -40,109 +41,66 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2 || nargout > 2) + if (args.length () != 2 || nargout > 2) { print_usage ("hess"); return retval; } - tree_constant arg = args(1).make_numeric (); + tree_constant arg = args(1); - int a_nr = arg.rows (); - int a_nc = arg.columns (); + int nr = arg.rows (); + int nc = arg.columns (); - if (a_nr == 0 || a_nc == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - warning ("hess: argument is empty matrix"); - Matrix m; - retval.resize (2); - retval(0) = m; - retval(1) = m; - } - else - error ("hess: empty matrix is invalid as argument"); + if (empty_arg ("hess", nr, nc) < 0) + return retval; - return retval; - } - - if (a_nr != a_nc) + if (nr != nc) { gripe_square_matrix_required ("hess"); return retval; } - Matrix tmp; - ComplexMatrix ctmp; + if (arg.is_real_type ()) + { + Matrix tmp = arg.matrix_value (); - if (arg.is_real_matrix ()) - { - tmp = arg.matrix_value (); - - HESS result (tmp); + if (! error_state) + { + HESS result (tmp); - if (nargout == 0 || nargout == 1) - { - retval.resize (1); - retval(0) = result.hess_matrix (); - } - else - { - retval.resize (2); - retval(0) = result.unitary_hess_matrix (); - retval(1) = result.hess_matrix (); + if (nargout == 0 || nargout == 1) + { + retval.resize (1); + retval(0) = result.hess_matrix (); + } + else + { + retval.resize (2); + retval(0) = result.unitary_hess_matrix (); + retval(1) = result.hess_matrix (); + } } } - else if (arg.is_complex_matrix ()) + else if (arg.is_complex_type ()) { - ctmp = arg.complex_matrix_value (); - ComplexHESS result (ctmp); + ComplexMatrix ctmp = arg.complex_matrix_value (); - if (nargout == 0 || nargout == 1) - { - retval.resize (1); - retval(0) = result.hess_matrix (); - } - else - { - retval.resize (2); - retval(0) = result.unitary_hess_matrix (); - retval(1) = result.hess_matrix (); - } - } - else if (arg.is_real_scalar ()) - { - double d = arg.double_value (); - if (nargout == 0 || nargout == 1) + if (! error_state) { - retval.resize (1); - retval(0) = d; - } - else - { - retval.resize (2); - retval(0) = 1; - retval(1) = d; - } - } - else if (arg.is_complex_scalar ()) - { - Complex c = arg.complex_value (); - if (nargout == 0 || nargout == 1) - { - retval.resize (1); - retval(0) = c; - } - else - { - retval.resize (2); - retval(0) = 1; - retval(1) = c; + ComplexHESS result (ctmp); + + if (nargout == 0 || nargout == 1) + { + retval.resize (1); + retval(0) = result.hess_matrix (); + } + else + { + retval.resize (2); + retval(0) = result.unitary_hess_matrix (); + retval(1) = result.hess_matrix (); + } } } else
--- a/src/ifft.cc +++ b/src/ifft.cc @@ -32,6 +32,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -40,51 +41,40 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2) + if (args.length () != 2) { print_usage ("ifft"); return retval; } - tree_constant tmp = args(1).make_numeric (); + tree_constant arg = args(1); - if (tmp.rows () == 0 || tmp.columns () == 0) + if (empty_arg ("ifft", arg.rows (), arg.columns ()) < 0) + return retval; + + if (arg.is_real_type ()) { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) + Matrix m = arg.matrix_value (); + + if (! error_state) { - if (flag < 0) - gripe_empty_arg ("ifft", 0); - - retval.resize (1, Matrix ()); + ComplexMatrix mifft = m.ifourier (); + retval = mifft; } - else - gripe_empty_arg ("ifft", 1); - - return retval; } - - if (tmp.is_real_matrix ()) + else if (arg.is_complex_type ()) { - Matrix m = tmp.matrix_value (); - ComplexMatrix mifft = m.ifourier (); - retval = mifft; - } - else if (tmp.is_complex_matrix ()) - { - ComplexMatrix m = tmp.complex_matrix_value (); - ComplexMatrix mifft = m.ifourier (); - retval = mifft; - } - else if (tmp.is_scalar_type ()) - { - error ("ifft: invalid scalar arguement"); + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix mifft = m.ifourier (); + retval = mifft; + } } else { - gripe_wrong_type_arg ("ifft", tmp); + gripe_wrong_type_arg ("ifft", arg); } return retval;
--- a/src/input.cc +++ b/src/input.cc @@ -753,13 +753,13 @@ char *prompt = "debug> "; if (nargin > 1) { - if (args(1).is_string ()) - prompt = args(1).string_value (); - else - { - error ("input: unrecognized argument"); - return retval; - } + prompt = args(1).string_value (); + + if (error_state) + { + error ("input: unrecognized argument"); + return retval; + } } again:
--- a/src/inv.cc +++ b/src/inv.cc @@ -32,6 +32,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -48,70 +49,59 @@ return retval; } - tree_constant tmp = args(1).make_numeric (); + tree_constant arg = args(1); + + int nr = arg.rows (); + int nc = arg.columns (); - int nr = tmp.rows (); - int nc = tmp.columns (); - if (nr == 0 || nc == 0) + if (empty_arg ("inverse", nr, nc) < 0) + return retval; + + if (nr != nc) { - int flag = user_pref.propagate_empty_matrices; - if (flag < 0) - gripe_empty_arg ("inverse", 0); - else if (flag == 0) - gripe_empty_arg ("inverse", 1); + gripe_square_matrix_required ("inverse"); + return retval; } - Matrix mtmp; - if (nr == 0 && nc == 0) - return mtmp; + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); - if (tmp.is_real_matrix ()) - { - Matrix m = tmp.matrix_value (); - if (m.rows () == m.columns ()) + if (! error_state) { int info; double rcond = 0.0; + Matrix minv = m.inverse (info, rcond); + if (info == -1) warning ("inverse: matrix singular to machine precision,\ rcond = %g", rcond); else retval = minv; } - else - gripe_square_matrix_required ("inverse"); } - else if (tmp.is_real_scalar ()) + else if (arg.is_complex_type ()) { - double d = 1.0 / tmp.double_value (); - retval = d; - } - else if (tmp.is_complex_matrix ()) - { - ComplexMatrix m = tmp.complex_matrix_value (); - if (m.rows () == m.columns ()) + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) { int info; double rcond = 0.0; + ComplexMatrix minv = m.inverse (info, rcond); + if (info == -1) warning ("inverse: matrix singular to machine precision,\ rcond = %g", rcond); else retval = minv; } - else - gripe_square_matrix_required ("inverse"); - } - else if (tmp.is_complex_scalar ()) - { - Complex c = 1.0 / tmp.complex_value (); - retval = c; } else { - gripe_wrong_type_arg ("inv", tmp); + gripe_wrong_type_arg ("inv", arg); } return retval;
--- a/src/log.cc +++ b/src/log.cc @@ -31,6 +31,7 @@ #include "user-prefs.h" #include "error.h" #include "gripes.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -41,94 +42,20 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2) + if (args.length () != 2) { print_usage ("logm"); return retval; } - tree_constant tmp = args(1).make_numeric ();; + tree_constant arg = args(1); - if (tmp.rows () == 0 || tmp.columns () == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - gripe_empty_arg ("logm", 0); - - retval.resize (1, Matrix ()); - return retval; - } - else - gripe_empty_arg ("logm", 1); - } - - if (tmp.is_real_matrix ()) - { - Matrix m = tmp.matrix_value (); - - int nr = m.rows (); - int nc = m.columns (); - - if (nr == 0 || nc == 0 || nr != nc) - gripe_square_matrix_required ("logm"); - else - { - EIG m_eig (m); - ComplexColumnVector lambda (m_eig.eigenvalues ()); - ComplexMatrix Q (m_eig.eigenvectors ()); + if (empty_arg ("logm", arg.rows (), arg.columns ()) < 0) + return retval; - for (int i = 0; i < nr; i++) - { - Complex elt = lambda.elem (i); - if (imag (elt) == 0.0 && real (elt) > 0.0) - lambda.elem (i) = log (real (elt)); - else - lambda.elem (i) = log (elt); - } - - ComplexDiagMatrix D (lambda); - ComplexMatrix result = Q * D * Q.inverse (); - - retval(0) = result; - } - } - else if (tmp.is_complex_matrix ()) + if (arg.is_real_scalar ()) { - ComplexMatrix m = tmp.complex_matrix_value (); - - int nr = m.rows (); - int nc = m.columns (); - - if (nr == 0 || nc == 0 || nr != nc) - gripe_square_matrix_required ("logm"); - else - { - EIG m_eig (m); - ComplexColumnVector lambda (m_eig.eigenvalues ()); - ComplexMatrix Q (m_eig.eigenvectors ()); - - for (int i = 0; i < nr; i++) - { - Complex elt = lambda.elem (i); - if (imag (elt) == 0.0 && real (elt) > 0.0) - lambda.elem (i) = log (real (elt)); - else - lambda.elem (i) = log (elt); - } - - ComplexDiagMatrix D (lambda); - ComplexMatrix result = Q * D * Q.inverse (); - - retval(0) = result; - } - } - else if (tmp.is_real_scalar ()) - { - double d = tmp.double_value (); + double d = arg.double_value (); if (d > 0.0) retval(0) = log (d); else @@ -137,14 +64,80 @@ retval(0) = log (dtmp); } } - else if (tmp.is_complex_scalar ()) + else if (arg.is_complex_scalar ()) + { + Complex c = arg.complex_value (); + retval(0) = log (c); + } + else if (arg.is_real_type ()) { - Complex c = tmp.complex_value (); - retval(0) = log (c); + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("logm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = log (real (elt)); + else + lambda.elem (i) = log (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("logm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = log (real (elt)); + else + lambda.elem (i) = log (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } } else { - gripe_wrong_type_arg ("logm", tmp); + gripe_wrong_type_arg ("logm", arg); } return retval; @@ -155,94 +148,20 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2) + if (args.length () != 2) { print_usage ("sqrtm"); return retval; } - tree_constant tmp = args(1).make_numeric ();; + tree_constant arg = args(1); - if (tmp.rows () == 0 || tmp.columns () == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - gripe_empty_arg ("sqrtm", 0); - - retval.resize (1, Matrix ()); - return retval; - } - else - gripe_empty_arg ("sqrtm", 1); - } - - if (tmp.is_real_matrix ()) - { - Matrix m = tmp.matrix_value (); - - int nr = m.rows (); - int nc = m.columns (); - - if (nr == 0 || nc == 0 || nr != nc) - gripe_square_matrix_required ("sqrtm"); - else - { - EIG m_eig (m); - ComplexColumnVector lambda (m_eig.eigenvalues ()); - ComplexMatrix Q (m_eig.eigenvectors ()); + if (empty_arg ("sqrtm", arg.rows (), arg.columns ())) + return retval; - for (int i = 0; i < nr; i++) - { - Complex elt = lambda.elem (i); - if (imag (elt) == 0.0 && real (elt) > 0.0) - lambda.elem (i) = sqrt (real (elt)); - else - lambda.elem (i) = sqrt (elt); - } - - ComplexDiagMatrix D (lambda); - ComplexMatrix result = Q * D * Q.inverse (); - - retval(0) = result; - } - } - else if (tmp.is_complex_matrix ()) + if (arg.is_real_scalar ()) { - ComplexMatrix m = tmp.complex_matrix_value (); - - int nr = m.rows (); - int nc = m.columns (); - - if (nr == 0 || nc == 0 || nr != nc) - gripe_square_matrix_required ("sqrtm"); - else - { - EIG m_eig (m); - ComplexColumnVector lambda (m_eig.eigenvalues ()); - ComplexMatrix Q (m_eig.eigenvectors ()); - - for (int i = 0; i < nr; i++) - { - Complex elt = lambda.elem (i); - if (imag (elt) == 0.0 && real (elt) > 0.0) - lambda.elem (i) = sqrt (real (elt)); - else - lambda.elem (i) = sqrt (elt); - } - - ComplexDiagMatrix D (lambda); - ComplexMatrix result = Q * D * Q.inverse (); - - retval(0) = result; - } - } - else if (tmp.is_real_scalar ()) - { - double d = tmp.double_value (); + double d = arg.double_value (); if (d > 0.0) retval(0) = sqrt (d); else @@ -251,14 +170,80 @@ retval(0) = sqrt (dtmp); } } - else if (tmp.is_complex_scalar ()) + else if (arg.is_complex_scalar ()) + { + Complex c = arg.complex_value (); + retval(0) = log (c); + } + else if (arg.is_real_type ()) { - Complex c = tmp.complex_value (); - retval(0) = log (c); + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("sqrtm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = sqrt (real (elt)); + else + lambda.elem (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("sqrtm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = sqrt (real (elt)); + else + lambda.elem (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } } else { - gripe_wrong_type_arg ("sqrtm", tmp); + gripe_wrong_type_arg ("sqrtm", arg); } return retval;
--- a/src/lsode.cc +++ b/src/lsode.cc @@ -84,7 +84,7 @@ { retval = tmp(0).vector_value (); - if (retval.length () == 0) + if (error_state || retval.length () == 0) gripe_user_supplied_eval ("lsode"); } else @@ -119,12 +119,34 @@ return retval; ColumnVector state = args(2).vector_value (); + + if (error_state) + { + error ("lsode: expecting state vector as second argument"); + return retval; + } + ColumnVector out_times = args(3).vector_value (); + + if (error_state) + { + error ("lsode: expecting output time vector as third argument"); + return retval; + } + ColumnVector crit_times; + int crit_times_set = 0; if (nargin > 4) { crit_times = args(4).vector_value (); + + if (error_state) + { + error ("lsode: expecting critical time vector as fourth argument"); + return retval; + } + crit_times_set = 1; } @@ -267,20 +289,25 @@ if (nargin == 1) { print_lsode_option_list (); + return retval; } else if (nargin == 3) { - if (args(1).is_string ()) + char *keyword = args(1).string_value (); + + if (! error_state) { - char *keyword = args(1).string_value (); double val = args(2).double_value (); - do_lsode_option (keyword, val); + + if (! error_state) + { + do_lsode_option (keyword, val); + return retval; + } } - else - print_usage ("lsode_options"); } - else - print_usage ("lsode_options"); + + print_usage ("lsode_options"); return retval; }
--- a/src/lu.cc +++ b/src/lu.cc @@ -31,6 +31,7 @@ #include "tree-const.h" #include "user-prefs.h" #include "gripes.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -39,37 +40,34 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2 || nargout > 3) + if (args.length () != 2 || nargout > 3) { print_usage ("lu"); return retval; } - tree_constant tmp = args(1).make_numeric ();; - - if (tmp.rows () == 0 || tmp.columns () == 0) + tree_constant arg = args(1); + + int nr = arg.rows (); + int nc = arg.columns (); + + if (empty_arg ("lu", nr, nc) < 0) + return retval; + + if (nr != nc) { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - gripe_empty_arg ("lu", 0); - - retval.resize (3, Matrix ()); - return retval; - } - else - gripe_empty_arg ("lu", 1); + gripe_square_matrix_required ("lu"); + return retval; } - if (tmp.is_real_matrix ()) + if (arg.is_real_type ()) { - Matrix m = tmp.matrix_value (); - if (m.rows () == m.columns ()) + Matrix m = arg.matrix_value (); + + if (! error_state) { LU fact (m); + switch (nargout) { case 0: @@ -82,6 +80,7 @@ retval(0) = L; } break; + case 3: default: retval(2) = fact.P (); @@ -90,15 +89,15 @@ break; } } - else - gripe_square_matrix_required ("lu"); } - else if (tmp.is_complex_matrix ()) + else if (arg.is_complex_matrix ()) { - ComplexMatrix m = tmp.complex_matrix_value (); - if (m.rows () == m.columns ()) + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) { ComplexLU fact (m); + switch (nargout) { case 0: @@ -111,6 +110,7 @@ retval(0) = L; } break; + case 3: default: retval(2) = fact.P (); @@ -119,26 +119,10 @@ break; } } - else - gripe_square_matrix_required ("lu"); - } - else if (tmp.is_real_scalar ()) - { - double d = tmp.double_value (); - retval(2) = 1.0; - retval(1) = d; - retval(0) = 1.0; - } - else if (tmp.is_complex_scalar ()) - { - Complex c = tmp.complex_value (); - retval(2) = 1.0; - retval(1) = c; - retval(0) = 1.0; } else { - gripe_wrong_type_arg ("lu", tmp); + gripe_wrong_type_arg ("lu", arg); } return retval;
--- a/src/minmax.cc +++ b/src/minmax.cc @@ -162,10 +162,10 @@ switch (nargin) { case 3: - arg2 = args(2).make_numeric (); + arg2 = args(2); // Fall through... case 2: - arg1 = args(1).make_numeric (); + arg1 = args(1); break; default: panic_impossible (); @@ -182,21 +182,29 @@ { retval(0) = arg1.complex_value (); } - else if (arg1.is_real_matrix ()) + else if (arg1.is_real_type ()) { Matrix m = arg1.matrix_value (); - if (m.rows () == 1) - retval(0) = m.row_min (); - else - retval(0) = tree_constant (m.column_min (), 0); + + if (! error_state) + { + if (m.rows () == 1) + retval(0) = m.row_min (); + else + retval(0) = tree_constant (m.column_min (), 0); + } } - else if (arg1.is_complex_matrix ()) + else if (arg1.is_complex_type ()) { ComplexMatrix m = arg1.complex_matrix_value (); - if (m.rows () == 1) - retval(0) = m.row_min (); - else - retval(0) = tree_constant (m.column_min (), 0); + + if (! error_state) + { + if (m.rows () == 1) + retval(0) = m.row_min (); + else + retval(0) = tree_constant (m.column_min (), 0); + } } else { @@ -216,32 +224,40 @@ retval(1) = 1; retval(0) = arg1.complex_value (); } - else if (arg1.is_real_matrix ()) + else if (arg1.is_real_type ()) { Matrix m = arg1.matrix_value (); - if (m.rows () == 1) + + if (! error_state) { - retval(1) = m.row_min_loc (); - retval(0) = m.row_min (); - } - else - { - retval(1) = tree_constant (m.column_min_loc (), 0); - retval(0) = tree_constant (m.column_min (), 0); + if (m.rows () == 1) + { + retval(1) = m.row_min_loc (); + retval(0) = m.row_min (); + } + else + { + retval(1) = tree_constant (m.column_min_loc (), 0); + retval(0) = tree_constant (m.column_min (), 0); + } } } - else if (arg1.is_complex_matrix ()) + else if (arg1.is_complex_type ()) { ComplexMatrix m = arg1.complex_matrix_value (); - if (m.rows () == 1) + + if (! error_state) { - retval(1) = m.row_min_loc (); - retval(0) = m.row_min (); - } - else - { - retval(1) = tree_constant (m.column_min_loc (), 0); - retval(0) = tree_constant (m.column_min (), 0); + if (m.rows () == 1) + { + retval(1) = m.row_min_loc (); + retval(0) = m.row_min (); + } + else + { + retval(1) = tree_constant (m.column_min_loc (), 0); + retval(0) = tree_constant (m.column_min (), 0); + } } } else @@ -255,38 +271,37 @@ if (arg1.rows () == arg2.rows () && arg1.columns () == arg2.columns ()) { -// XXX FIXME XXX -- I don't think this is quite right. - if (arg1.is_real_scalar ()) - { - double result; - double a_elem = arg1.double_value (); - double b_elem = arg2.double_value (); - result = MIN (a_elem, b_elem); - retval(0) = result; - } - else if (arg1.is_complex_scalar ()) + if (arg1.is_real_type () && arg2.is_real_type ()) { - Complex result; - Complex a_elem = arg1.complex_value (); - Complex b_elem = arg2.complex_value (); - if (abs (a_elem) < abs (b_elem)) - result = a_elem; - else - result = b_elem; - retval(0) = result; + Matrix m1 = arg1.matrix_value (); + + if (! error_state) + { + Matrix m2 = arg2.matrix_value (); + + if (! error_state) + { + Matrix result = min (m1, m2); + if (! error_state) + retval(0) = result; + } + } } - else if (arg1.is_real_matrix ()) + else if (arg1.is_complex_matrix () || arg2.is_complex_type ()) { - Matrix result; - result = min (arg1.matrix_value (), arg2.matrix_value ()); - retval(0) = result; - } - else if (arg1.is_complex_matrix ()) - { - ComplexMatrix result; - result = min (arg1.complex_matrix_value (), - arg2.complex_matrix_value ()); - retval(0) = result; + ComplexMatrix m1 = arg1.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix m2 = arg2.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix result = min (m1, m2); + if (! error_state) + retval(0) = result; + } + } } else { @@ -322,10 +337,10 @@ switch (nargin) { case 3: - arg2 = args(2).make_numeric (); + arg2 = args(2); // Fall through... case 2: - arg1 = args(1).make_numeric (); + arg1 = args(1); break; default: panic_impossible ();
--- a/src/npsol.cc +++ b/src/npsol.cc @@ -159,7 +159,7 @@ { retval = tmp(0).vector_value (); - if (retval.length () <= 0) + if (error_state || retval.length () <= 0) error ("npsol: error evaluating constraints"); } else @@ -296,7 +296,7 @@ ColumnVector x = args(1).vector_value (); - if (x.capacity () == 0) + if (error_state || x.capacity () == 0) { error ("npsol: expecting vector as first argument"); return retval; @@ -319,7 +319,8 @@ int lb_len = lb.capacity (); int ub_len = ub.capacity (); - if (lb_len != ub_len || lb_len != x.capacity ()) + + if (error_state || lb_len != ub_len || lb_len != x.capacity ()) { error ("npsol: lower and upper bounds and decision variable vector"); error ("must all have the same number of elements"); @@ -366,15 +367,22 @@ if (! npsol_constraints) { ColumnVector lub = args(nargin-1).vector_value (); - Matrix c = args(nargin-2).matrix_value (); ColumnVector llb = args(nargin-3).vector_value (); - if (llb.capacity () == 0 || lub.capacity () == 0) + if (error_state || llb.capacity () == 0 || lub.capacity () == 0) { error ("npsol: bounds for linear constraints must be vectors"); return retval; } + Matrix c = args(nargin-2).matrix_value (); + + if (error_state) + { + error ("npsol: invalid linear constraint matrix"); + return retval; + } + if (! linear_constraints_ok (x, llb, c, lub, "npsol", 1)) return retval; @@ -405,12 +413,12 @@ ColumnVector nlub = args(nargin-1).vector_value (); ColumnVector nllb = args(nargin-3).vector_value (); - NLFunc const_func (npsol_constraint_function); - - if (! nonlinear_constraints_ok - (x, nllb, npsol_constraint_function, nlub, "npsol", 1)) + if (error_state + || (! nonlinear_constraints_ok + (x, nllb, npsol_constraint_function, nlub, "npsol", 1))) return retval; + NLFunc const_func (npsol_constraint_function); NLConst nonlinear_constraints (nllb, const_func, nlub); if (nargin == 6) @@ -448,23 +456,30 @@ ColumnVector nlub = args(nargin-1).vector_value (); ColumnVector nllb = args(nargin-3).vector_value (); - NLFunc const_func (npsol_constraint_function); - - if (! nonlinear_constraints_ok - (x, nllb, npsol_constraint_function, nlub, "npsol", 1)) + if (error_state + || (! nonlinear_constraints_ok + (x, nllb, npsol_constraint_function, nlub, "npsol", 1))) return retval; + NLFunc const_func (npsol_constraint_function); NLConst nonlinear_constraints (nllb, const_func, nlub); ColumnVector lub = args(nargin-4).vector_value (); - Matrix c = args(nargin-5).matrix_value (); ColumnVector llb = args(nargin-6).vector_value (); - if (llb.capacity () == 0 || lub.capacity () == 0) + if (error_state || llb.capacity () == 0 || lub.capacity () == 0) { error ("npsol: bounds for linear constraints must be vectors"); return retval; } + + Matrix c = args(nargin-5).matrix_value (); + + if (error_state) + { + error ("npsol: invalid linear constraint matrix"); + return retval; + } if (! linear_constraints_ok (x, llb, c, lub, "npsol", 1)) return retval; @@ -750,22 +765,25 @@ if (nargin == 1) { print_npsol_option_list (); + return retval; } else if (nargin == 3) { - if (args(1).is_string ()) + char *keyword = args(1).string_value (); + + if (! error_state) { - char *keyword = args(1).string_value (); double val = args(2).double_value (); - do_npsol_option (keyword, val); + + if (! error_state) + { + do_npsol_option (keyword, val); + return retval; + } } - else - print_usage ("npsol_options"); } - else - { - print_usage ("npsol_options"); - } + + print_usage ("npsol_options"); #endif
--- a/src/octave.cc +++ b/src/octave.cc @@ -34,7 +34,6 @@ #include <sys/stat.h> #include <time.h> #include <pwd.h> -#include <setjmp.h> #include <stdlib.h> #include <string.h> #include <signal.h> @@ -43,6 +42,11 @@ #include <strstream.h> #include <fstream.h> +extern "C" +{ +#include <setjmp.h> +} + #include "getopt.h" #include "lo-error.h" @@ -469,6 +473,8 @@ initialize_pager (); + install_signal_handlers (); + initialize_history (); initialize_file_io (); @@ -479,8 +485,6 @@ initialize_readline (); - install_signal_handlers (); - if (! inhibit_startup_message) cout << "Octave, version " << version_string << ". Copyright (C) 1992, 1993, 1994 John W. Eaton.\n" @@ -738,14 +742,14 @@ tree_constant eval_string (const tree_constant& arg, int& parse_status) { - if (! arg.is_string ()) + char *string = arg.string_value (); + + if (error_state) { error ("eval: expecting string argument"); return -1; } - char *string = arg.string_value (); - // Yes Virginia, we always print here... return eval_string (string, 1, 1, parse_status); @@ -786,9 +790,15 @@ tree_constant tc_command = args(1); - if (tc_command.is_string ()) + char *tmp_str = tc_command.string_value (); + + if (error_state) { - iprocstream cmd (tc_command.string_value ()); + error ("shell_cmd: expecting string as first argument"); + } + else + { + iprocstream cmd (tmp_str); ostrstream output_buf; @@ -812,8 +822,6 @@ else maybe_page_output (output_buf); } - else - error ("shell_cmd: expecting string as first argument"); return retval; }
--- a/src/pt-const.h +++ b/src/pt-const.h @@ -204,16 +204,14 @@ tree_constant all (void) const { return rep->all (); } tree_constant any (void) const { return rep->any (); } -// Broader classifications. - - int is_scalar_type (void) const { return rep->is_scalar_type (); } - int is_matrix_type (void) const { return rep->is_matrix_type (); } - int is_real_type (void) const { return rep->is_real_type (); } int is_complex_type (void) const { return rep->is_complex_type (); } -// These need better names, since a range really is a numeric type. +// Would be nice to get rid of the next four functions: + + int is_scalar_type (void) const { return rep->is_scalar_type (); } + int is_matrix_type (void) const { return rep->is_matrix_type (); } int is_numeric_type (void) const { return rep->is_numeric_type (); } @@ -284,43 +282,6 @@ void convert_to_row_or_column_vector (void) { rep->convert_to_row_or_column_vector (); } -// These need better names, since a range really is a numeric type. - - void force_numeric (int force_str_conv = 0) - { rep->force_numeric (force_str_conv); } - - tree_constant make_numeric (int force_str_conv = 0) const - { - if (is_numeric_type ()) - return *this; - else - return rep->make_numeric (force_str_conv); - } - - tree_constant make_numeric_or_range (void) const - { - if (is_numeric_type () || is_range ()) - return *this; - else - return rep->make_numeric (); - } - - tree_constant make_numeric_or_magic (void) const - { - if (is_numeric_type () || is_magic_colon ()) - return *this; - else - return rep->make_numeric (); - } - - tree_constant make_numeric_or_range_or_magic (void) const - { - if (is_numeric_type () || is_range () || is_magic_colon ()) - return *this; - else - return rep->make_numeric (); - } - // Increment or decrement this constant. void bump_value (tree_expression::type et) @@ -402,6 +363,49 @@ { return rep->const_type (); } // ------------------------------------------------------------------- + +private: + +// Can we make these go away? + +// These need better names, since a range really is a numeric type. + + void force_numeric (int force_str_conv = 0) + { rep->force_numeric (force_str_conv); } + + tree_constant make_numeric (int force_str_conv = 0) const + { + if (is_numeric_type ()) + return *this; + else + return rep->make_numeric (force_str_conv); + } + +#if 0 + tree_constant make_numeric_or_range (void) const + { + if (is_numeric_type () || is_range ()) + return *this; + else + return rep->make_numeric (); + } +#endif + + tree_constant make_numeric_or_magic (void) const + { + if (is_numeric_type () || is_magic_colon ()) + return *this; + else + return rep->make_numeric (); + } + + tree_constant make_numeric_or_range_or_magic (void) const + { + if (is_numeric_type () || is_range () || is_magic_colon ()) + return *this; + else + return rep->make_numeric (); + } }; // XXX FIXME XXX -- this is not used very much now. Perhaps it can be
--- a/src/pt-exp-base.cc +++ b/src/pt-exp-base.cc @@ -314,8 +314,7 @@ goto done; } - Octave_object otmp = elem->eval (0); - tree_constant tmp = otmp(0); + tree_constant tmp = elem->eval (0); if (error_state || tmp.is_undefined ()) { retval = tree_constant (); @@ -496,17 +495,7 @@ str_ptr += nc; } } - else if (tmp.is_range ()) - { - tmp.force_numeric (1); - if (tmp.is_real_scalar ()) - m (put_row, put_col) = tmp.double_value (); - else if (tmp.is_real_matrix ()) - m.insert (tmp.matrix_value (), put_row, put_col); - else - panic_impossible (); - } - else if (tmp.is_real_matrix ()) + else if (tmp.is_real_matrix () || tmp.is_range ()) { cm.insert (tmp.matrix_value (), put_row, put_col); } @@ -537,17 +526,7 @@ str_ptr += nc; } } - else if (tmp.is_range ()) - { - tmp.force_numeric (1); - if (tmp.is_real_scalar ()) - m (put_row, put_col) = tmp.double_value (); - else if (tmp.is_real_matrix ()) - m.insert (tmp.matrix_value (), put_row, put_col); - else - panic_impossible (); - } - else if (tmp.is_real_matrix ()) + else if (tmp.is_real_matrix () || tmp.is_range ()) { m.insert (tmp.matrix_value (), put_row, put_col); } @@ -1234,8 +1213,7 @@ case tree_expression::transpose: if (op) { - Octave_object tmp = op->eval (0); - tree_constant u = tmp(0); + tree_constant u = op->eval (0); if (error_state) eval_error (); else if (u.is_defined ()) @@ -1347,14 +1325,12 @@ case tree_expression::or: if (op1) { - Octave_object tmp = op1->eval (0); - tree_constant a = tmp(0); + tree_constant a = op1->eval (0); if (error_state) eval_error (); else if (a.is_defined () && op2) { - tmp = op2->eval (0); - tree_constant b = tmp (0); + tree_constant b = op2->eval (0); if (error_state) eval_error (); else if (b.is_defined ()) @@ -1376,8 +1352,7 @@ int result = 0; if (op1) { - Octave_object tmp = op1->eval (0); - tree_constant a = tmp(0); + tree_constant a = op1->eval (0); if (error_state) { eval_error (); @@ -1410,8 +1385,7 @@ if (op2) { - tmp = op2->eval (0); - tree_constant b = tmp(0); + tree_constant b = op2->eval (0); if (error_state) { eval_error (); @@ -1535,8 +1509,7 @@ if (rhs) { - Octave_object tmp = rhs->eval (0); - tree_constant rhs_val = tmp(0); + tree_constant rhs_val = rhs->eval (0); if (error_state) { if (error_state) @@ -1844,8 +1817,7 @@ if (error_state || ! op1 || ! op2) return retval; - Octave_object otmp = op1->eval (0); - tree_constant tmp = otmp(0); + tree_constant tmp = op1->eval (0); if (tmp.is_undefined ()) { @@ -1853,16 +1825,15 @@ return retval; } - tmp = tmp.make_numeric (); - if (! tmp.is_scalar_type ()) + double base = tmp.double_value (); + + if (error_state) { - eval_error ("base for colon expression must be a scalar"); + eval_error ("evaluating colon expression"); return retval; } - double base = tmp.double_value (); - - otmp = op2->eval (0); - tmp = otmp(0); + + tmp = op2->eval (0); if (tmp.is_undefined ()) { @@ -1870,19 +1841,18 @@ return retval; } - tmp = tmp.make_numeric (); - if (! tmp.is_scalar_type ()) + double limit = tmp.double_value (); + + if (error_state) { - eval_error ("limit for colon expression must be a scalar"); + eval_error ("evaluating colon expression"); return retval; } - double limit = tmp.double_value (); double inc = 1.0; if (op3) { - otmp = op3->eval (0); - tmp = otmp(0); + tmp = op3->eval (0); if (tmp.is_undefined ()) { @@ -1890,14 +1860,13 @@ return retval; } - tmp = tmp.make_numeric (); - if (! tmp.is_scalar_type ()) + inc = tmp.double_value (); + + if (error_state) { - eval_error ("increment for colon expression must be a scalar"); + eval_error ("evaluating colon expression"); return retval; } - else - inc = tmp.double_value (); } retval = tree_constant (base, limit, inc);
--- a/src/qpsol.cc +++ b/src/qpsol.cc @@ -97,21 +97,24 @@ } ColumnVector x = args(1).vector_value (); - if (x.capacity () == 0) + + if (error_state || x.capacity () == 0) { error ("qpsol: expecting vector as first argument"); return retval; } Matrix H = args(2).matrix_value (); - if (H.rows () != H.columns () || H.rows () != x.capacity ()) + + if (error_state || H.rows () != H.columns () || H.rows () != x.capacity ()) { error ("qpsol: H must be a square matrix consistent with the size of x"); return retval; } ColumnVector c = args(3).vector_value (); - if (c.capacity () != x.capacity ()) + + if (error_state || c.capacity () != x.capacity ()) { error ("qpsol: c must be a vector the same size as x"); return retval; @@ -125,7 +128,8 @@ int lb_len = lb.capacity (); int ub_len = ub.capacity (); - if (lb_len != ub_len || lb_len != x.capacity ()) + + if (error_state || lb_len != ub_len || lb_len != x.capacity ()) { error ("qpsol: lower and upper bounds and decision variable vector"); error ("must all have the same number of elements"); @@ -167,15 +171,22 @@ if (nargin == 7 || nargin == 9) { ColumnVector lub = args(nargin-1).vector_value (); - Matrix A = args(nargin-2).matrix_value (); ColumnVector llb = args(nargin-3).vector_value (); - if (llb.capacity () == 0 || lub.capacity () == 0) + if (error_state || llb.capacity () == 0 || lub.capacity () == 0) { error ("qpsol: bounds for linear constraints must be vectors"); return retval; } + Matrix A = args(nargin-2).matrix_value (); + + if (error_state) + { + error ("qpsol: invalid linear constraint matrix"); + return retval; + } + if (! linear_constraints_ok (x, llb, A, lub, "qpsol", 1)) return retval; @@ -360,22 +371,25 @@ if (nargin == 1) { print_qpsol_option_list (); + return retval; } else if (nargin == 3) { - if (args(1).is_string ()) + char *keyword = args(1).string_value (); + + if (! error_state) { - char *keyword = args(1).string_value (); double val = args(2).double_value (); - do_qpsol_option (keyword, val); + + if (! error_state) + { + do_qpsol_option (keyword, val); + return retval; + } } - else - print_usage ("qpsol_options"); } - else - { - print_usage ("qpsol_options"); - } + + print_usage ("qpsol_options"); #endif
--- a/src/qr.cc +++ b/src/qr.cc @@ -34,6 +34,7 @@ #include "tree-const.h" #include "user-prefs.h" #include "gripes.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -65,90 +66,59 @@ return retval; } - tree_constant tmp = args(1).make_numeric (); + tree_constant arg = args(1); - if (tmp.rows () == 0 || tmp.columns () == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - gripe_empty_arg ("qr", 0); - Matrix m; - retval(2) = m; - retval(1) = m; - retval(0) = m; - } - else - gripe_empty_arg ("qr", 1); - - return retval; - } + if (empty_arg ("qr", arg.rows (), arg.columns ()) < 0) + return retval; QR::type type = nargout == 1 ? QR::raw : (nargin == 3 ? QR::economy : QR::std); - if (tmp.is_real_matrix ()) + if (arg.is_real_type ()) { - Matrix m = tmp.matrix_value (); - if (nargout < 3) + Matrix m = arg.matrix_value (); + + if (! error_state) { - QR fact (m, type); - retval(1) = fact.R (); - retval(0) = fact.Q (); - } - else - { - QRP fact (m, type); - retval(2) = fact.P (); - retval(1) = fact.R (); - retval(0) = fact.Q (); + if (nargout < 3) + { + QR fact (m, type); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + else + { + QRP fact (m, type); + retval(2) = fact.P (); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } } } - else if (tmp.is_complex_matrix ()) - { - ComplexMatrix m = tmp.complex_matrix_value (); - if (nargout < 3) - { - ComplexQR fact (m, type); - retval(1) = fact.R (); - retval(0) = fact.Q (); - } - else - { - ComplexQRP fact (m, type); - retval(2) = fact.P (); - retval(1) = fact.R (); - retval(0) = fact.Q (); - } - } - else if (tmp.is_real_scalar ()) + else if (arg.is_complex_type ()) { - double d = tmp.double_value (); - if (nargout == 1) - retval(0) = d; - else + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) { - retval(2) = 1.0; - retval(1) = d; - retval(0) = 1.0; - } - } - else if (tmp.is_complex_scalar ()) - { - Complex c = tmp.complex_value (); - if (nargout == 1) - retval(0) = c; - else - { - retval(2) = 1.0; - retval(1) = c; - retval(0) = 1.0; + if (nargout < 3) + { + ComplexQR fact (m, type); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + else + { + ComplexQRP fact (m, type); + retval(2) = fact.P (); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } } } else { - gripe_wrong_type_arg ("qr", tmp); + gripe_wrong_type_arg ("qr", arg); } return retval;
--- a/src/quad.cc +++ b/src/quad.cc @@ -66,7 +66,15 @@ } if (tmp.length () && tmp(0).is_defined ()) - retval = tmp(0).double_value (); + { + retval = tmp(0).double_value (); + + if (error_state) + { + quad_integration_error = 1; // XXX FIXME XXX + gripe_user_supplied_eval ("quad"); + } + } else { quad_integration_error = 1; // XXX FIXME XXX @@ -107,8 +115,21 @@ return retval; double a = args(2).double_value (); + + if (error_state) + { + error ("quad: expecting second argument to be a scalar"); + return retval; + } + double b = args(3).double_value (); + if (error_state) + { + error ("quad: expecting third argument to be a scalar"); + return retval; + } + int indefinite = 0; IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite; double bound = 0.0; @@ -147,21 +168,40 @@ error("quad: singularities not allowed on infinite intervals"); return retval; } + have_sing = 1; + sing = args(5).vector_value (); + + if (error_state) + { + error ("quad: expecting vector of singularities as fourth argument"); + return retval; + } + case 5: tol = args(4).vector_value (); + + if (error_state) + { + error ("quad: expecting vector of tolerances as fifth argument"); + return retval; + } + switch (tol.capacity ()) { case 2: reltol = tol.elem (1); + case 1: abstol = tol.elem (0); break; + default: error ("quad: expecting tol to contain no more than two values"); return retval; } + case 4: if (indefinite) { @@ -185,17 +225,16 @@ } } break; + default: panic_impossible (); break; } - retval.resize (4); - + retval(3) = abserr; + retval(2) = nfun; + retval(1) = ier; retval(0) = val; - retval(1) = ier; - retval(2) = nfun; - retval(3) = abserr; return retval; } @@ -299,20 +338,27 @@ int nargin = args.length (); if (nargin == 1) - print_quad_option_list (); + { + print_quad_option_list (); + return retval; + } else if (nargin == 3) { - if (args(1).is_string ()) + char *keyword = args(1).string_value (); + + if (! error_state) { - char *keyword = args(1).string_value (); double val = args(2).double_value (); - do_quad_option (keyword, val); + + if (! error_state) + { + do_quad_option (keyword, val); + return retval; + } } - else - print_usage ("quad_options"); } - else - print_usage ("quad_options"); + + print_usage ("quad_options"); return retval; }
--- a/src/qzval.cc +++ b/src/qzval.cc @@ -38,6 +38,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -61,107 +62,115 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 3 || nargout > 1) + if (args.length () != 3 || nargout > 1) { print_usage ("qzvalue"); return retval; } - tree_constant arga = args(1).make_numeric (); - tree_constant argb = args(2).make_numeric(); + tree_constant arg_a = args(1); + tree_constant arg_b = args(2); + + int a_nr = arg_a.rows(); + int a_nc = arg_a.columns(); - if (arga.is_empty () || argb.is_empty ()) - retval = vector_of_empties (nargout, "qzvalue"); - else - { + int b_nr = arg_b.rows(); + int b_nc = arg_b.columns(); + + if (empty_arg ("qzvalue", a_nr, a_nc) < 0 + || empty_arg ("qzvalue", b_nr, b_nc) < 0) + return retval; // Arguments are not empty, so check for correct dimensions. - int a_rows = arga.rows(); - int a_cols = arga.columns(); - int b_rows = argb.rows(); - int b_cols = argb.columns(); - - if ((a_rows != a_cols) || (b_rows != b_cols)) - { - gripe_square_matrix_required ("qzvalue: first two parameters:"); - return retval; - } - else if (a_rows != b_rows) - { - gripe_nonconformant (); - return retval; - } + if (a_nr != a_nc || b_nr != b_nc) + { + gripe_square_matrix_required ("qzvalue: first two parameters:"); + return retval; + } + + if (a_nr != b_nr) + { + gripe_nonconformant (); + return retval; + } // Dimensions look o.k., let's solve the problem. - retval.resize (nargout ? nargout : 1); - - if (arga.is_complex_type () || argb.is_complex_type ()) - error ("qzvalue: cannot yet do complex matrix arguments\n"); - else - { + if (arg_a.is_complex_type () || arg_b.is_complex_type ()) + { + error ("qzvalue: cannot yet do complex matrix arguments\n"); + return retval; + } // Do everything in real arithmetic. - Matrix jnk (a_rows, a_rows, 0.0); + Matrix jnk (a_nr, a_nr, 0.0); - ColumnVector alfr (a_rows); - ColumnVector alfi (a_rows); - ColumnVector beta (a_rows); + ColumnVector alfr (a_nr); + ColumnVector alfi (a_nr); + ColumnVector beta (a_nr); - long matz = 0; - int info; + long matz = 0; + int info; // XXX FIXME ??? XXX - double eps = DBL_EPSILON; + double eps = DBL_EPSILON; + + Matrix ca = arg_a.matrix_value (); - Matrix ca = arga.matrix_value (); - Matrix cb = argb.matrix_value (); + if (error_state) + return retval; + + Matrix cb = arg_b.matrix_value (); + + if (error_state) + return retval; // Use EISPACK qz functions. - F77_FCN (qzhes) (&a_rows, &a_rows, ca.fortran_vec (), - cb.fortran_vec (), &matz, jnk.fortran_vec ()); + F77_FCN (qzhes) (&a_nr, &a_nr, ca.fortran_vec (), + cb.fortran_vec (), &matz, jnk.fortran_vec ()); - F77_FCN (qzit) (&a_rows, &a_rows, ca.fortran_vec (), - cb.fortran_vec (), &eps, &matz, - jnk.fortran_vec (), &info); + F77_FCN (qzit) (&a_nr, &a_nr, ca.fortran_vec (), + cb.fortran_vec (), &eps, &matz, + jnk.fortran_vec (), &info); - if (info) - error ("qzvalue: trouble in qzit, info = %d", info); + if (info) + error ("qzvalue: trouble in qzit, info = %d", info); - F77_FCN (qzval) (&a_rows, &a_rows, ca.fortran_vec (), - cb.fortran_vec (), alfr.fortran_vec (), - alfi.fortran_vec (), beta.fortran_vec (), - &matz, jnk.fortran_vec ()); + F77_FCN (qzval) (&a_nr, &a_nr, ca.fortran_vec (), + cb.fortran_vec (), alfr.fortran_vec (), + alfi.fortran_vec (), beta.fortran_vec (), + &matz, jnk.fortran_vec ()); // Count and extract finite generalized eigenvalues. - int i, cnt; - Complex Im (0, 1); - for (i = 0, cnt = 0; i < a_rows; i++) - if (beta (i) != 0) - cnt++; + int i; + int cnt = 0; + + Complex Im (0, 1); - ComplexColumnVector cx (cnt, 0.0); + for (i = 0; i < a_nr; i++) + if (beta (i) != 0) + cnt++; - for (i = 0; i < a_rows; i++) - { - if (beta (i) != 0) - { + ComplexColumnVector cx (cnt, 0.0); + + for (i = 0; i < a_nr; i++) + { + if (beta (i) != 0) + { // Finite generalized eigenvalue. - cnt--; - cx (cnt) = (alfr (i) + Im * alfi (i)) / beta (i); - } - } - retval(0) = cx; + cnt--; + cx (cnt) = (alfr (i) + Im * alfi (i)) / beta (i); } } + + retval = cx; + return retval; }
--- a/src/rand.cc +++ b/src/rand.cc @@ -177,9 +177,10 @@ } else if (tmp.is_scalar_type ()) { - n = NINT (tmp.double_value ()); - m = n; - goto gen_matrix; + m = n = NINT (tmp.double_value ()); + + if (! error_state) + goto gen_matrix; } else if (tmp.is_range ()) { @@ -206,13 +207,21 @@ && strcmp (args(1).string_value (), "seed") == 0) { double d = args(2).double_value (); - set_rand_seed (d); + + if (! error_state) + set_rand_seed (d); } else { n = NINT (args(1).double_value ()); - m = NINT (args(2).double_value ()); - goto gen_matrix; + + if (! error_state) + { + m = NINT (args(2).double_value ()); + + if (! error_state) + goto gen_matrix; + } } }
--- a/src/schur.cc +++ b/src/schur.cc @@ -30,6 +30,7 @@ #include "tree-const.h" #include "user-prefs.h" +#include "utils.h" #include "error.h" #include "gripes.h" #include "help.h" @@ -58,111 +59,77 @@ return retval; } - tree_constant arg = args(1).make_numeric (); + tree_constant arg = args(1); + + char *ord = "U"; + if (nargin == 3) + { + ord = args(2).string_value (); - char *ord; - if (nargin != 3) - ord = "U"; - else - ord = args(2).string_value (); + if (error_state) + { + error ("schur: expecting string as third argument"); + return retval; + } + } if (*ord != 'U' && *ord != 'A' && *ord != 'D' && *ord != 'u' && *ord != 'a' && *ord != 'd') { warning ("schur: incorrect ordered schur argument `%c'", *ord); - Matrix m; - retval.resize (2); - retval(0) = m; - retval(1) = m; return retval; } - int a_nr = arg.rows (); - int a_nc = arg.columns (); + + int nr = arg.rows (); + int nc = arg.columns (); - if (a_nr == 0 || a_nc == 0) - { - int flag = user_pref.propagate_empty_matrices; - if (flag != 0) - { - if (flag < 0) - warning ("schur: argument is empty matrix"); - Matrix m; - retval.resize (2); - retval(0) = m; - retval(1) = m; - } - else - error ("schur: empty matrix is invalid as argument"); + if (empty_arg ("schur", nr, nc) < 0) + return retval; - return retval; - } - if (a_nr != a_nc) + if (nr != nc) { gripe_square_matrix_required ("schur"); return retval; } - Matrix tmp; - ComplexMatrix ctmp; - if (arg.is_real_matrix ()) { - tmp = arg.matrix_value (); + Matrix tmp = arg.matrix_value (); - SCHUR result (tmp,ord); - - if (nargout == 0 || nargout == 1) + if (! error_state) { - retval(0) = result.schur_matrix (); - } - else - { - retval(1) = result.schur_matrix (); - retval(0) = result.unitary_matrix (); + SCHUR result (tmp,ord); + + if (nargout == 0 || nargout == 1) + { + retval(0) = result.schur_matrix (); + } + else + { + retval(1) = result.schur_matrix (); + retval(0) = result.unitary_matrix (); + } } } else if (arg.is_complex_matrix ()) { - ctmp = arg.complex_matrix_value (); + ComplexMatrix ctmp = arg.complex_matrix_value (); - ComplexSCHUR result (ctmp,ord); - - if (nargout == 0 || nargout == 1) + if (! error_state) { - retval(0) = result.schur_matrix (); - } - else - { - retval(1) = result.schur_matrix (); - retval(0) = result.unitary_matrix (); + ComplexSCHUR result (ctmp,ord); + + if (nargout == 0 || nargout == 1) + { + retval(0) = result.schur_matrix (); + } + else + { + retval(1) = result.schur_matrix (); + retval(0) = result.unitary_matrix (); + } } } - else if (arg.is_real_scalar ()) - { - double d = arg.double_value (); - if (nargout == 0 || nargout == 1) - { - retval(0) = d; - } - else - { - retval(1) = d; - retval(0) = 1.0; - } - } - else if (arg.is_complex_scalar ()) - { - Complex c = arg.complex_value (); - if (nargout == 0 || nargout == 1) - { - retval(0) = c; - } - else - { - retval(1) = c; - retval(0) = 1.0; - } - } else { gripe_wrong_type_arg ("schur", arg);
--- a/src/sort.cc +++ b/src/sort.cc @@ -184,77 +184,73 @@ else retval.resize (1); - tree_constant tmp = args(1); + tree_constant arg = args(1); - if (tmp.is_real_scalar ()) + if (arg.is_real_type ()) { - retval(0) = tmp.double_value (); - if (return_idx) - retval(1) = 1.0; - } - else if (tmp.is_complex_scalar ()) - { - retval(0) = tmp.complex_value (); - if (return_idx) - retval(1) = 1.0; - } - else if (tmp.is_real_matrix () || tmp.is_string () || tmp.is_range ()) - { - Matrix m = tmp.matrix_value (); - if (m.rows () == 1) + Matrix m = arg.matrix_value (); + + if (! error_state) { - int nc = m.columns (); - RowVector v (nc); - for (int i = 0; i < nc; i++) - v.elem (i) = m.elem (0, i); - RowVector idx; - mx_sort (v, idx, return_idx); + if (m.rows () == 1) + { + int nc = m.columns (); + RowVector v (nc); + for (int i = 0; i < nc; i++) + v.elem (i) = m.elem (0, i); + RowVector idx; + mx_sort (v, idx, return_idx); - retval(0) = tree_constant (v, 0); - if (return_idx) - retval(1) = tree_constant (idx, 0); - } - else - { + retval(0) = tree_constant (v, 0); + if (return_idx) + retval(1) = tree_constant (idx, 0); + } + else + { // Sorts m in place, optionally computes index Matrix. - Matrix idx; - mx_sort (m, idx, return_idx); + Matrix idx; + mx_sort (m, idx, return_idx); - retval(0) = m; - if (return_idx) - retval(1) = idx; + retval(0) = m; + if (return_idx) + retval(1) = idx; + } } } - else if (tmp.is_complex_matrix ()) + else if (arg.is_complex_type ()) { - ComplexMatrix cm = tmp.complex_matrix_value (); - if (cm.rows () == 1) + ComplexMatrix cm = arg.complex_matrix_value (); + + if (! error_state) { - int nc = cm.columns (); - ComplexRowVector cv (nc); - for (int i = 0; i < nc; i++) - cv.elem (i) = cm.elem (0, i); - RowVector idx; - mx_sort (cv, idx, return_idx); + if (cm.rows () == 1) + { + int nc = cm.columns (); + ComplexRowVector cv (nc); + for (int i = 0; i < nc; i++) + cv.elem (i) = cm.elem (0, i); + RowVector idx; + mx_sort (cv, idx, return_idx); - retval(0) = tree_constant (cv, 0); - if (return_idx) - retval(1) = tree_constant (idx, 0); - } - else - { + retval(0) = tree_constant (cv, 0); + if (return_idx) + retval(1) = tree_constant (idx, 0); + } + else + { // Sorts cm in place, optionally computes index Matrix. - Matrix idx; - mx_sort (cm, idx, return_idx); + Matrix idx; + mx_sort (cm, idx, return_idx); - retval(0) = cm; - if (return_idx) - retval(1) = idx; + retval(0) = cm; + if (return_idx) + retval(1) = idx; + } } } else { - gripe_wrong_type_arg ("sort", tmp); + gripe_wrong_type_arg ("sort", arg); } return retval;
--- a/src/svd.cc +++ b/src/svd.cc @@ -43,7 +43,7 @@ argument, an `economy' sized factorization is computed that omits\n\ unnecessary rows and columns of U and V") { - Octave_object retval (3, Matrix ()); + Octave_object retval; int nargin = args.length (); @@ -53,52 +53,55 @@ return retval; } - SVD::type type = (nargin == 3) ? SVD::economy : SVD::std; + tree_constant arg = args(1); - tree_constant arg = args(1); + if (empty_arg ("svd", arg.rows (), arg.columns ()) < 0) + return retval; + + SVD::type type = (nargin == 3) ? SVD::economy : SVD::std; if (arg.is_real_type ()) { Matrix tmp = arg.matrix_value (); - if (error_state || empty_arg ("svd", tmp.rows (), tmp.columns ())) - return retval; + if (! error_state) + { + SVD result (tmp, type); - SVD result (tmp, type); - - DiagMatrix sigma = result.singular_values (); + DiagMatrix sigma = result.singular_values (); - if (nargout == 0 || nargout == 1) - { - retval(0) = tree_constant (sigma.diag (), 1); - } - else - { - retval(2) = result.right_singular_matrix (); - retval(1) = sigma; - retval(0) = result.left_singular_matrix (); + if (nargout == 0 || nargout == 1) + { + retval(0) = tree_constant (sigma.diag (), 1); + } + else + { + retval(2) = result.right_singular_matrix (); + retval(1) = sigma; + retval(0) = result.left_singular_matrix (); + } } } else if (arg.is_complex_type ()) { ComplexMatrix ctmp = arg.complex_matrix_value (); - if (error_state || empty_arg ("svd", ctmp.rows (), ctmp.columns ())) - return retval; + if (! error_state) + { + ComplexSVD result (ctmp, type); - ComplexSVD result (ctmp, type); - - DiagMatrix sigma = result.singular_values (); + DiagMatrix sigma = result.singular_values (); - if (nargout == 0 || nargout == 1) - { - retval(0) = tree_constant (sigma.diag (), 1); - } - else - { - retval(2) = result.right_singular_matrix (); - retval(1) = sigma; - retval(0) = result.left_singular_matrix (); + if (nargout == 0 || nargout == 1) + { + retval(0) = tree_constant (sigma.diag (), 1); + } + else + { + retval(2) = result.right_singular_matrix (); + retval(1) = sigma; + retval(0) = result.left_singular_matrix (); + } } } else
--- a/src/syl.cc +++ b/src/syl.cc @@ -37,6 +37,7 @@ #include "user-prefs.h" #include "gripes.h" #include "error.h" +#include "utils.h" #include "help.h" #include "defun-dld.h" @@ -60,57 +61,67 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 4 || nargout > 1) + if (args.length () != 4 || nargout > 1) { print_usage ("syl"); return retval; } - tree_constant arga = args(1).make_numeric (); - tree_constant argb = args(2).make_numeric (); - tree_constant argc = args(3).make_numeric (); + tree_constant arg_a = args(1); + tree_constant arg_b = args(2); + tree_constant arg_c = args(3); + + int a_nr = arg_a.rows (); + int a_nc = arg_a.columns (); - if (arga.is_empty () || argb.is_empty () || argc.is_empty ()) - retval = vector_of_empties (nargout, "syl"); - else - { + int b_nr = arg_b.rows (); + int b_nc = arg_b.columns (); + + int c_nr = arg_c.rows (); + int c_nc = arg_c.columns (); + + if (empty_arg ("syl", a_nr, a_nc) < 0 + || empty_arg ("syl", b_nr, b_nc) < 0 + || empty_arg ("syl", c_nr, c_nc) < 0) + return retval; // Arguments are not empty, so check for correct dimensions. - int a_rows = arga.rows (); - int a_cols = arga.columns (); - int b_rows = argb.rows (); - int b_cols = argb.columns (); - int c_rows = argc.rows (); - int c_cols = argc.columns (); - - if ((a_rows != a_cols) || (b_rows != b_cols)) - { - gripe_square_matrix_required ("syl: first two parameters:"); - return retval; - } - else if ((a_rows != c_rows) || (b_rows != c_cols)) - { - gripe_nonconformant (); - return retval; - } + if (a_nr != a_nc || b_nr != b_nc) + { + gripe_square_matrix_required ("syl: first two parameters:"); + return retval; + } + else if (a_nr != c_nr || b_nr != c_nc) + { + gripe_nonconformant (); + return retval; + } // Dimensions look o.k., let's solve the problem. - retval.resize (nargout ? nargout : 1); - - if (arga.is_complex_type () || argb.is_complex_type () - || argc.is_complex_type ()) + if (arg_a.is_complex_type () + || arg_b.is_complex_type () + || arg_c.is_complex_type ()) { // Do everything in complex arithmetic; - ComplexMatrix ca = arga.complex_matrix_value (); - ComplexMatrix cb = argb.complex_matrix_value (); - ComplexMatrix cc = argc.complex_matrix_value (); - + ComplexMatrix ca = arg_a.complex_matrix_value (); + + if (error_state) + return retval; + + ComplexMatrix cb = arg_b.complex_matrix_value (); + + if (error_state) + return retval; + + ComplexMatrix cc = arg_c.complex_matrix_value (); + + if (error_state) + return retval; + // Compute Schur decompositions ComplexSCHUR as (ca, "U"); @@ -131,25 +142,36 @@ int info; int one = 1; - F77_FCN (ztrsyl) ("N", "N", &one, &a_rows, &b_rows, - sch_a.fortran_vec (), &a_rows, - sch_b.fortran_vec (), &b_rows, - cx.fortran_vec (), &a_rows, &scale, &info, + F77_FCN (ztrsyl) ("N", "N", &one, &a_nr, &b_nr, + sch_a.fortran_vec (), &a_nr, + sch_b.fortran_vec (), &b_nr, + cx.fortran_vec (), &a_nr, &scale, &info, 1L, 1L); cx = -ua * cx * ub.hermitian (); - retval(0) = cx; + retval = cx; } else { // Do everything in real arithmetic; - Matrix ca = arga.matrix_value (); - Matrix cb = argb.matrix_value (); - Matrix cc = argc.matrix_value (); - + Matrix ca = arg_a.matrix_value (); + + if (error_state) + return retval; + + Matrix cb = arg_b.matrix_value (); + + if (error_state) + return retval; + + Matrix cc = arg_c.matrix_value (); + + if (error_state) + return retval; + // Compute Schur decompositions. SCHUR as (ca, "U"); @@ -170,10 +192,10 @@ int info; int one = 1; - F77_FCN (dtrsyl) ("N", "N", &one, &a_rows, &b_rows, - sch_a.fortran_vec (), &a_rows, - sch_b.fortran_vec (), &b_rows, - cx.fortran_vec (), &a_rows, &scale, &info, + F77_FCN (dtrsyl) ("N", "N", &one, &a_nr, &b_nr, + sch_a.fortran_vec (), &a_nr, + sch_b.fortran_vec (), &b_nr, + cx.fortran_vec (), &a_nr, &scale, &info, 1L, 1L); if (info) @@ -181,9 +203,9 @@ cx = -ua*cx*ub.transpose (); - retval(0) = cx; + retval = cx; } - } + return retval; }
--- a/src/sysdep.cc +++ b/src/sysdep.cc @@ -396,13 +396,18 @@ int nargin = args.length (); - if (nargin == 2 && args(1).is_string ()) + if (nargin == 2) { - char *value = getenv (args(1).string_value ()); - if (value) - retval = value; - else - retval = ""; + char *name = args(1).string_value (); + + if (! error_state) + { + char *value = getenv (name); + if (value) + retval = value; + else + retval = ""; + } } else print_usage ("getenv"); @@ -448,13 +453,17 @@ { case 2: { - int delay = NINT (args(1).double_value ()); - if (delay > 0) + double dval = args(1).double_value (); + + if (! error_state) { - sleep (delay); - break; + int delay = NINT (dval); + if (delay > 0) + sleep (delay); } } + break; + default: if (kbhit () == EOF) clean_up_and_exit (0);
--- a/src/tc-rep.cc +++ b/src/tc-rep.cc @@ -904,8 +904,6 @@ retval = matrix->elem (0, 0); else gripe_invalid_conversion ("real matrix", "real scalar"); - - retval = octave_NaN; } break; @@ -1053,7 +1051,7 @@ Complex TC_REP::complex_value (int force_string_conversion) const { - Complex retval; + Complex retval (octave_NaN, octave_NaN); switch (type_tag) { @@ -1077,8 +1075,6 @@ } else gripe_invalid_conversion ("real matrix", "real scalar"); - - retval = octave_NaN; } break; @@ -1185,8 +1181,13 @@ char * TC_REP::string_value (void) const { - assert (type_tag == string_constant); - return string; + if (type_tag == string_constant) + return string; + else + { + gripe_invalid_conversion (type_as_string (), "string"); + return 0; + } } Range
--- a/src/tc-rep.h +++ b/src/tc-rep.h @@ -118,18 +118,6 @@ tree_constant all (void) const; tree_constant any (void) const; - int is_scalar_type (void) const - { - return (type_tag == scalar_constant - || type_tag == complex_scalar_constant); - } - - int is_matrix_type (void) const - { - return (type_tag == matrix_constant - || type_tag == complex_matrix_constant); - } - int is_real_type (void) const { return (type_tag == scalar_constant @@ -144,6 +132,20 @@ || type_tag == complex_scalar_constant); } +// Would be nice to get rid of the next four functions: + + int is_scalar_type (void) const + { + return (type_tag == scalar_constant + || type_tag == complex_scalar_constant); + } + + int is_matrix_type (void) const + { + return (type_tag == matrix_constant + || type_tag == complex_matrix_constant); + } + int is_numeric_type (void) const { return (type_tag == scalar_constant
--- a/src/utils.cc +++ b/src/utils.cc @@ -45,7 +45,6 @@ #include <unistd.h> #endif #include <sys/param.h> -#include <setjmp.h> #include <string.h> #include <limits.h> #include <iostream.h> @@ -54,12 +53,14 @@ #include <Complex.h> -#ifndef HAVE_STRNCASECMP extern "C" { +#include <setjmp.h> + +#ifndef HAVE_STRNCASECMP extern int strncasecmp (const char*, const char*, size_t); +#endif } -#endif extern "C" { @@ -514,14 +515,18 @@ if (nr == 0 || nc == 0) { - is_empty = 0; - int flag = user_pref.propagate_empty_matrices; if (flag < 0) - gripe_empty_arg (name, 0); - else if (flag > 0) - gripe_empty_arg (name, 1); + { + gripe_empty_arg (name, 0); + is_empty = 1; + } + else if (is_empty > 0) + { + gripe_empty_arg (name, 1); + is_empty = -1; + } } return is_empty;
--- a/src/variables.cc +++ b/src/variables.cc @@ -133,15 +133,15 @@ { tree_fvc *ans = 0; - if (! arg.is_string ()) + char *fcn_name = arg.string_value (); + + if (error_state) { if (warn) error ("%s: expecting function name as argument", warn_for); return ans; } - char *fcn_name = arg.string_value (); - symbol_record *sr = 0; if (fcn_name) sr = lookup_by_name (fcn_name); @@ -184,9 +184,7 @@ { Octave_object retval = 0.0; - int nargin = args.length (); - - if (nargin != 2 || ! args(1).is_string ()) + if (args.length () != 2) { print_usage ("is_global"); return retval; @@ -194,6 +192,12 @@ char *name = args(1).string_value (); + if (error_state) + { + error ("is_global: expecting string argument"); + return retval; + } + symbol_record *sr = curr_sym_tab->lookup (name, 0, 0); retval = (double) (sr && sr->is_linked_to_global ()); @@ -209,9 +213,7 @@ { Octave_object retval; - int nargin = args.length (); - - if (nargin != 2 || ! args(1).is_string ()) + if (args.length () != 2) { print_usage ("exist"); return retval; @@ -219,6 +221,12 @@ char *name = args(1).string_value (); + if (error_state) + { + error ("exist: expecting string argument"); + return retval; + } + symbol_record *sr = curr_sym_tab->lookup (name, 0, 0); if (! sr) sr = global_sym_tab->lookup (name, 0, 0); @@ -678,6 +686,7 @@ if (! error_state && val.is_string ()) { char *s = val.string_value (); + if (s) retval = strsave (s); }