Mercurial > hg > octave-nkf
changeset 1357:749071f48336
[project @ 1995-09-05 20:28:59 by jwe]
author | jwe |
---|---|
date | Tue, 05 Sep 1995 20:28:59 +0000 |
parents | 8b51c1738882 |
children | dc9c01f66a19 |
files | src/balance.cc src/eig.cc src/expm.cc src/find.cc src/givens.cc src/log.cc src/lpsolve.cc src/minmax.cc src/npsol.cc src/qpsol.cc src/qzval.cc src/rand.cc src/sort.cc src/syl.cc |
diffstat | 14 files changed, 104 insertions(+), 88 deletions(-) [+] |
line wrap: on
line diff
--- a/src/balance.cc +++ b/src/balance.cc @@ -75,8 +75,8 @@ 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 ()) { @@ -94,7 +94,7 @@ int a_nr = arg_a.rows (); int a_nc = arg_a.columns (); -// Check argument 1 dimensions. + // Check argument 1 dimensions. int arg_is_empty = empty_arg ("balance", a_nr, a_nc); @@ -109,7 +109,7 @@ return retval; } -// Extract argument 1 parameter for both AEP and GEP. + // Extract argument 1 parameter for both AEP and GEP. Matrix aa; ComplexMatrix caa; @@ -121,13 +121,13 @@ if (error_state) return retval; -// Treat AEP/ GEP cases. + // Treat AEP/GEP cases. switch (my_nargin) { case 1: -// Algebraic eigenvalue problem. + // Algebraic eigenvalue problem. if (arg_a.is_complex_type ()) { @@ -157,16 +157,16 @@ case 2: { -// Generalized eigenvalue problem. + // Generalized eigenvalue problem. -// 1st we have to check argument 2 dimensions and type... + // 1st we have to check argument 2 dimensions and type... tree_constant arg_b = args(1); int b_nr = arg_b.rows (); int b_nc = arg_b.columns (); -// Check argument 2 dimensions -- must match arg 1. + // Check argument 2 dimensions -- must match arg 1. if (b_nr != b_nc || b_nr != a_nr) { @@ -174,8 +174,8 @@ return retval; } -// Now, extract the second matrix... -// Extract argument 1 parameter for both AEP and GEP. + // Now, extract the second matrix... + // Extract argument 1 parameter for both AEP and GEP. Matrix bb; ComplexMatrix cbb; @@ -187,7 +187,7 @@ if (error_state) return retval; -// Both matrices loaded, now let's check what kind of arithmetic: + // Both matrices loaded, now let's check what kind of arithmetic: if (arg_a.is_complex_type () || arg_b.is_complex_type ()) { @@ -197,8 +197,8 @@ if (arg_b.is_real_type ()) cbb = bb; -// Compute magnitudes of elements for balancing purposes. -// Surely there's a function I can call someplace! + // 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_nc; j++)
--- a/src/eig.cc +++ b/src/eig.cc @@ -99,7 +99,7 @@ } else { -// Blame it on Matlab. + // Blame it on Matlab. ComplexDiagMatrix d (result.eigenvalues ());
--- a/src/expm.cc +++ b/src/expm.cc @@ -70,7 +70,7 @@ tree_constant arg = args(0); -// Constants for matrix exponential calculation. + // Constants for matrix exponential calculation. static double padec [] = { @@ -110,8 +110,7 @@ if (arg.is_real_type ()) { - -// Compute the exponential. + // Compute the exponential. Matrix m = arg.matrix_value (); @@ -120,7 +119,7 @@ double trshift = 0; // trace shift value -// Preconditioning step 1: trace normalization. + // Preconditioning step 1: trace normalization. for (i = 0; i < nc; i++) trshift += m.elem (i, i); @@ -128,13 +127,13 @@ for (i = 0; i < nc; i++) m.elem (i, i) -= trshift; -// Preconditioning step 2: balancing. + // Preconditioning step 2: balancing. AEPBALANCE mbal (m, balance_job); m = mbal.balanced_matrix (); Matrix d = mbal.balancing_matrix (); -// Preconditioning step 3: scaling. + // Preconditioning step 3: scaling. ColumnVector work(nc); inf_norm = F77_FCN (dlange, DLANGE) ("I", nc, nc, @@ -143,7 +142,7 @@ sqpow = (int) (1.0 + log (inf_norm) / log (2.0)); -// Check whether we need to square at all. + // Check whether we need to square at all. if (sqpow < 0) sqpow = 0; @@ -155,12 +154,12 @@ m = m / inf_norm; } -// npp, dpp: pade' approx polynomial matrices. + // npp, dpp: pade' approx polynomial matrices. Matrix npp (nc, nc, 0.0); Matrix dpp = npp; -// now powers a^8 ... a^1. + // Now powers a^8 ... a^1. minus_one_j = -1; for (j = 7; j >= 0; j--) @@ -169,7 +168,8 @@ dpp = m * dpp + m * (minus_one_j * padec[j]); minus_one_j *= -1; } -// Zero power. + + // Zero power. dpp = -dpp; for(j = 0; j < nc; j++) @@ -178,11 +178,11 @@ dpp.elem (j, j) += 1.0; } -// Compute pade approximation = inverse (dpp) * npp. + // Compute pade approximation = inverse (dpp) * npp. Matrix result = dpp.solve (npp); -// Reverse preconditioning step 3: repeated squaring. + // Reverse preconditioning step 3: repeated squaring. while (sqpow) { @@ -190,7 +190,7 @@ sqpow--; } -// Reverse preconditioning step 2: inverse balancing. + // Reverse preconditioning step 2: inverse balancing. result = result.transpose(); d = d.transpose (); @@ -198,7 +198,7 @@ result = d.solve (result); result = result.transpose (); -// Reverse preconditioning step 1: fix trace normalization. + // Reverse preconditioning step 1: fix trace normalization. result = result * exp (trshift); @@ -213,7 +213,7 @@ Complex trshift = 0.0; // trace shift value -// Preconditioning step 1: trace normalization. + // Preconditioning step 1: trace normalization. for (i = 0; i < nc; i++) trshift += m.elem (i, i); @@ -221,13 +221,13 @@ for (i = 0; i < nc; i++) m.elem (i, i) -= trshift; -// Preconditioning step 2: eigenvalue balancing. + // Preconditioning step 2: eigenvalue balancing. ComplexAEPBALANCE mbal (m, balance_job); m = mbal.balanced_matrix (); ComplexMatrix d = mbal.balancing_matrix (); -// Preconditioning step 3: scaling. + // Preconditioning step 3: scaling. ColumnVector work (nc); inf_norm = F77_FCN (zlange, ZLANGE) ("I", nc, nc, @@ -236,7 +236,7 @@ sqpow = (int) (1.0 + log (inf_norm) / log (2.0)); -// Check whether we need to square at all. + // Check whether we need to square at all. if (sqpow < 0) sqpow = 0; @@ -248,12 +248,12 @@ m = m / inf_norm; } -// npp, dpp: pade' approx polynomial matrices. + // npp, dpp: pade' approx polynomial matrices. ComplexMatrix npp (nc, nc, 0.0); ComplexMatrix dpp = npp; -// Now powers a^8 ... a^1. + // Now powers a^8 ... a^1. minus_one_j = -1; for (j = 7; j >= 0; j--) @@ -263,7 +263,7 @@ minus_one_j *= -1; } -// Zero power. + // Zero power. dpp = -dpp; for (j = 0; j < nc; j++) @@ -272,11 +272,11 @@ dpp.elem (j, j) += 1.0; } -// Compute pade approximation = inverse (dpp) * npp. + // Compute pade approximation = inverse (dpp) * npp. ComplexMatrix result = dpp.solve (npp); -// Reverse preconditioning step 3: repeated squaring. + // Reverse preconditioning step 3: repeated squaring. while (sqpow) { @@ -284,9 +284,9 @@ sqpow--; } -// reverse preconditioning step 2: inverse balancing XXX FIXME XXX: -// should probably do this with lapack calls instead of a complete -// matrix inversion. + // 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 (); @@ -294,7 +294,7 @@ result = d.solve (result); result = result.transpose (); -// Reverse preconditioning step 1: fix trace normalization. + // Reverse preconditioning step 1: fix trace normalization. result = result * exp (trshift);
--- a/src/find.cc +++ b/src/find.cc @@ -47,8 +47,8 @@ for (int i = 0; i < count; i++) tmp (i) = nr * (j_idx (i) - 1.0) + i_idx (i); -// If the original argument was a row vector, force a row vector of -// indices to be returned. + // If the original argument was a row vector, force a row + // vector of indices to be returned. retval(0) = tree_constant (tmp, (nr != 1)); } @@ -56,14 +56,18 @@ case 3: retval(2) = val; -// Fall through! + // 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)); + + // If you want this to work more like Matlab, use + // + // retval(0) = tree_constant (i_idx, (nr != 1)); + // + // instead of the previous statement. + break; default:
--- a/src/givens.cc +++ b/src/givens.cc @@ -91,7 +91,9 @@ if (error_state) return retval; - cx = x; // copy to complex just in case + // Convert to complex just in case... + + cx = x; } if (arg_b.is_complex_type ()) @@ -108,10 +110,12 @@ if (error_state) return retval; - cy = y; // copy to complex just in case + // Convert to complex just in case... + + cy = y; } -// Now compute the rotation. + // Now compute the rotation. double cc; if (arg_a.is_complex_type () || arg_b.is_complex_type ()) @@ -122,7 +126,7 @@ switch (nargout) { - case 0: // output a matrix + case 0: case 1: { ComplexMatrix g (2, 2); @@ -135,7 +139,7 @@ } break; - case 2: // output scalar values + case 2: retval(0) = cc; retval(1) = cs; break; @@ -153,7 +157,7 @@ switch (nargout) { - case 0: // output a matrix + case 0: case 1: { Matrix g (2, 2); @@ -166,7 +170,7 @@ } break; - case 2: // output scalar values + case 2: retval(0) = cc; retval(1) = s; break;
--- a/src/log.cc +++ b/src/log.cc @@ -35,7 +35,8 @@ #include "user-prefs.h" #include "utils.h" -// XXX FIXME XXX -- the next two functions should really be just one... +// XXX FIXME XXX -- the next two functions should really be just +// one... DEFUN_DLD_BUILTIN ("logm", Flogm, Slogm, 2, 1, "logm (X): matrix logarithm")
--- a/src/lpsolve.cc +++ b/src/lpsolve.cc @@ -37,7 +37,8 @@ { Octave_object retval; -// Force a bad value of inform, and empty matrices for x and phi. + // Force a bad value of inform, and empty matrices for x and phi. + Matrix m; retval(2) = -1.0; retval(1) = m;
--- a/src/minmax.cc +++ b/src/minmax.cc @@ -382,7 +382,7 @@ { case 2: arg2 = args(1); -// Fall through... + // Fall through... case 1: arg1 = args(0); @@ -610,7 +610,7 @@ { case 2: arg2 = args(1); -// Fall through... + // Fall through... case 1: arg1 = args(0);
--- a/src/npsol.cc +++ b/src/npsol.cc @@ -273,7 +273,8 @@ #if defined (NPSOL_MISSING) -// Force a bad value of inform, and empty matrices for x, phi, and lambda. + // Force a bad value of inform, and empty matrices for x, phi, and + // lambda. retval.resize (4, Matrix ()); @@ -445,6 +446,7 @@ if (! npsol_constraints) { // Produce error message. + is_valid_function (args(nargin-2), "npsol", 1); } else
--- a/src/qpsol.cc +++ b/src/qpsol.cc @@ -79,7 +79,8 @@ #if defined (QPSOL_MISSING) -// Force a bad value of inform, and empty matrices for x, phi, and lambda. + // Force a bad value of inform, and empty matrices for x, phi, and + // lambda. retval.resize (4, Matrix ());
--- a/src/qzval.cc +++ b/src/qzval.cc @@ -89,7 +89,7 @@ else if (arg_a_is_empty || arg_b_is_empty) return retval; -// Arguments are not empty, so check for correct dimensions. + // Arguments are not empty, so check for correct dimensions. if (a_nr != a_nc || b_nr != b_nc) { @@ -103,7 +103,7 @@ return retval; } -// Dimensions look o.k., let's solve the problem. + // Dimensions look o.k., let's solve the problem. if (arg_a.is_complex_type () || arg_b.is_complex_type ()) { @@ -111,7 +111,7 @@ return retval; } -// Do everything in real arithmetic. + // Do everything in real arithmetic. Matrix jnk (a_nr, a_nr, 0.0); @@ -122,7 +122,7 @@ long matz = 0; int info; -// XXX FIXME ??? XXX + // XXX FIXME ??? XXX double eps = DBL_EPSILON; Matrix ca = arg_a.matrix_value (); @@ -135,7 +135,7 @@ if (error_state) return retval; -// Use EISPACK qz functions. + // Use EISPACK qz functions. F77_FCN (qzhes, QZHES) (a_nr, a_nr, ca.fortran_vec (), cb.fortran_vec (), matz, @@ -153,7 +153,7 @@ alfi.fortran_vec (), beta.fortran_vec (), matz, jnk.fortran_vec ()); -// Count and extract finite generalized eigenvalues. + // Count and extract finite generalized eigenvalues. int i; int cnt = 0; @@ -170,8 +170,7 @@ { if (beta (i) != 0) { - -// Finite generalized eigenvalue. + // Finite generalized eigenvalue. cnt--; cx (cnt) = (alfr (i) + Im * alfi (i)) / beta (i);
--- a/src/rand.cc +++ b/src/rand.cc @@ -128,10 +128,11 @@ static int initialized = 0; if (! initialized) { -// Make the random number generator give us a different sequence every -// time we start octave unless we specifically set the seed. The -// technique used below will cycle monthly, but it it does seem to -// work ok to give fairly different seeds each time Octave starts. + // Make the random number generator give us a different sequence + // every time we start octave unless we specifically set the + // seed. The technique used below will cycle monthly, but it it + // does seem to work ok to give fairly different seeds each time + // Octave starts. #if 0 int s0 = 1234567890; @@ -214,7 +215,8 @@ } else if (tmp.is_matrix_type ()) { -// XXX FIXME XXX -- this should probably use the function from data.cc. + // XXX FIXME XXX -- this should probably use the function + // from data.cc. Matrix a = args(0).matrix_value ();
--- a/src/sort.cc +++ b/src/sort.cc @@ -209,7 +209,8 @@ } else { -// Sorts m in place, optionally computes index Matrix. + // Sorts m in place, optionally computes index Matrix. + Matrix idx; mx_sort (m, idx, return_idx); @@ -240,7 +241,8 @@ } else { -// Sorts cm in place, optionally computes index Matrix. + // Sorts cm in place, optionally computes index Matrix. + Matrix idx; mx_sort (cm, idx, return_idx);
--- a/src/syl.cc +++ b/src/syl.cc @@ -92,7 +92,7 @@ else if (arg_a_is_empty || arg_b_is_empty || arg_c_is_empty) return retval; -// Arguments are not empty, so check for correct dimensions. + // Arguments are not empty, so check for correct dimensions. if (a_nr != a_nc || b_nr != b_nc) { @@ -105,14 +105,13 @@ return retval; } -// Dimensions look o.k., let's solve the problem. + // Dimensions look o.k., let's solve the problem. if (arg_a.is_complex_type () || arg_b.is_complex_type () || arg_c.is_complex_type ()) { - -// Do everything in complex arithmetic; + // Do everything in complex arithmetic; ComplexMatrix ca = arg_a.complex_matrix_value (); @@ -129,12 +128,12 @@ if (error_state) return retval; -// Compute Schur decompositions + // Compute Schur decompositions ComplexSCHUR as (ca, "U"); ComplexSCHUR bs (cb, "U"); -// Transform cc to new coordinates. + // Transform cc to new coordinates. ComplexMatrix ua = as.unitary_matrix (); ComplexMatrix sch_a = as.schur_matrix (); @@ -143,7 +142,8 @@ ComplexMatrix cx = ua.hermitian () * cc * ub; -// Solve the sylvester equation, back-transform, and return the solution. + // Solve the sylvester equation, back-transform, and return + // the solution. double scale; int info; @@ -160,8 +160,7 @@ } else { - -// Do everything in real arithmetic; + // Do everything in real arithmetic. Matrix ca = arg_a.matrix_value (); @@ -178,12 +177,12 @@ if (error_state) return retval; -// Compute Schur decompositions. + // Compute Schur decompositions. SCHUR as (ca, "U"); SCHUR bs (cb, "U"); -// Transform cc to new coordinates. + // Transform cc to new coordinates. Matrix ua = as.unitary_matrix (); Matrix sch_a = as.schur_matrix (); @@ -192,7 +191,8 @@ Matrix cx = ua.transpose () * cc * ub; -// Solve the sylvester equation, back-transform, and return the solution. + // Solve the sylvester equation, back-transform, and return + // the solution. double scale; int info;