Author: dimpbx
Date: Wed Aug 11 15:55:14 2010
New Revision: 984453
URL: http://svn.apache.org/viewvc?rev=984453&view=rev
Log:
MathUtils.safeNorm (translation of enorm.f from Minpack) added
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java?rev=984453&r1=984452&r2=984453&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java
Wed Aug 11 15:55:14 2010
@@ -1955,4 +1955,119 @@ public static void checkOrder(double[] v
checkOrder(val, OrderDirection.DECREASING, strict);
}
}
+
+ /**
+ * Returns the Cartesian norm (2-norm), handling both overflow and
underflow.
+ * Translation of the minpack enorm subroutine.
+ *
+ * The redistribution policy for MINPACK is available <a
+ * href="http://www.netlib.org/minpack/disclaimer";>here</a>, for
convenience, it
+ * is reproduced below.</p>
+ *
+ * <table border="0" width="80%" cellpadding="10" align="center"
bgcolor="#E0E0E0">
+ * <tr><td>
+ * Minpack Copyright Notice (1999) University of Chicago.
+ * All rights reserved
+ * </td></tr>
+ * <tr><td>
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * <ol>
+ * <li>Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following
disclaimer.</li>
+ * <li>Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.</li>
+ * <li>The end-user documentation included with the redistribution, if
any,
+ * must include the following acknowledgment:
+ * <code>This product includes software developed by the University
of
+ * Chicago, as Operator of Argonne National Laboratory.</code>
+ * Alternately, this acknowledgment may appear in the software
itself,
+ * if and wherever such third-party acknowledgments normally
appear.</li>
+ * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
+ * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
+ * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
+ * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
+ * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
+ * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
+ * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
+ * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
+ * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
+ * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
+ * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
+ * BE CORRECTED.</strong></li>
+ * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
+ * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
+ * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
+ * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
+ * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
+ * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
+ * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
+ * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
+ * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
+ * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
+ * <ol></td></tr>
+ * </table>
+ *
+ * @param v vector of doubles
+ * @return the 2-norm of the vector
+ */
+ public static double safeNorm(double[] v) {
+ double rdwarf = 3.834e-20;
+ double rgiant = 1.304e+19;
+ double s1=0.0;
+ double s2=0.0;
+ double s3=0.0;
+ double x1max = 0.0;
+ double x3max = 0.0;
+ double floatn = (double)v.length;
+ double agiant = rgiant/floatn;
+ for (int i=0;i<v.length;i++) {
+ double xabs = Math.abs(v[i]);
+ if (xabs<rdwarf || xabs>agiant) {
+ if (xabs>rdwarf) {
+ if (xabs>x1max) {
+ double r=x1max/xabs;
+ s1=1.0+s1*r*r;
+ x1max=xabs;
+ } else {
+ double r=xabs/x1max;
+ s1+=r*r;
+ }
+ } else {
+ if (xabs>x3max) {
+ double r=x3max/xabs;
+ s3=1.0+s3*r*r;
+ x3max=xabs;
+ } else {
+ if (xabs!=0.0) {
+ double r=xabs/x3max;
+ s3+=r*r;
+ }
+ }
+ }
+ } else {
+ s2+=xabs*xabs;
+ }
+ }
+ double norm;
+ if (s1!=0.0) {
+ norm = x1max*Math.sqrt(s1+(s2/x1max)/x1max);
+ } else {
+ if (s2==0.0) {
+ norm = x3max*Math.sqrt(s3);
+ } else {
+ if (s2>=x3max) {
+ norm = Math.sqrt(s2*(1.0+(x3max/s2)*(x3max*s3)));
+ } else {
+ norm = Math.sqrt(x3max*((s2/x3max)+(x3max*s3)));
+ }
+ }
+ }
+ return norm;
+}
+
}
