[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [46820] branches/smoke2/intern/smoke/ intern: Forgot to commit rest of smoke
Daniel Genrich
daniel.genrich at gmx.net
Sun May 20 23:00:33 CEST 2012
Revision: 46820
http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=46820
Author: genscher
Date: 2012-05-20 21:00:33 +0000 (Sun, 20 May 2012)
Log Message:
-----------
Forgot to commit rest of smoke
Added Paths:
-----------
branches/smoke2/intern/smoke/intern/tnt/
branches/smoke2/intern/smoke/intern/tnt/jama_eig.h
branches/smoke2/intern/smoke/intern/tnt/jama_lu.h
branches/smoke2/intern/smoke/intern/tnt/tnt.h
branches/smoke2/intern/smoke/intern/tnt/tnt_array1d.h
branches/smoke2/intern/smoke/intern/tnt/tnt_array1d_utils.h
branches/smoke2/intern/smoke/intern/tnt/tnt_array2d.h
branches/smoke2/intern/smoke/intern/tnt/tnt_array2d_utils.h
branches/smoke2/intern/smoke/intern/tnt/tnt_array3d.h
branches/smoke2/intern/smoke/intern/tnt/tnt_array3d_utils.h
branches/smoke2/intern/smoke/intern/tnt/tnt_cmat.h
branches/smoke2/intern/smoke/intern/tnt/tnt_fortran_array1d.h
branches/smoke2/intern/smoke/intern/tnt/tnt_fortran_array1d_utils.h
branches/smoke2/intern/smoke/intern/tnt/tnt_fortran_array2d.h
branches/smoke2/intern/smoke/intern/tnt/tnt_fortran_array2d_utils.h
branches/smoke2/intern/smoke/intern/tnt/tnt_fortran_array3d.h
branches/smoke2/intern/smoke/intern/tnt/tnt_fortran_array3d_utils.h
branches/smoke2/intern/smoke/intern/tnt/tnt_i_refvec.h
branches/smoke2/intern/smoke/intern/tnt/tnt_math_utils.h
branches/smoke2/intern/smoke/intern/tnt/tnt_sparse_matrix_csr.h
branches/smoke2/intern/smoke/intern/tnt/tnt_stopwatch.h
branches/smoke2/intern/smoke/intern/tnt/tnt_subscript.h
branches/smoke2/intern/smoke/intern/tnt/tnt_vec.h
branches/smoke2/intern/smoke/intern/tnt/tnt_version.h
Added: branches/smoke2/intern/smoke/intern/tnt/jama_eig.h
===================================================================
--- branches/smoke2/intern/smoke/intern/tnt/jama_eig.h (rev 0)
+++ branches/smoke2/intern/smoke/intern/tnt/jama_eig.h 2012-05-20 21:00:33 UTC (rev 46820)
@@ -0,0 +1,1053 @@
+/** \file smoke/intern/tnt/jama_eig.h
+ * \ingroup smoke
+ */
+#ifndef JAMA_EIG_H
+#define JAMA_EIG_H
+
+
+#include "tnt_array1d.h"
+#include "tnt_array2d.h"
+#include "tnt_math_utils.h"
+
+#include <algorithm>
+// for min(), max() below
+
+#include <cmath>
+// for fabs() below
+
+using namespace TNT;
+using namespace std;
+
+// NT debugging
+//static int gEigenDebug=0;
+//if(gEigenDebug) std::cerr<<"n="<<n<<" m="<<m<<" l="<<l<<"\n";
+// m has to be smaller l! in line 262
+// gcc can get confused with abs calls, replaced by fabs
+
+namespace JAMA
+{
+
+/**
+
+ Computes eigenvalues and eigenvectors of a real (non-complex)
+ matrix.
+<P>
+ If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
+ diagonal and the eigenvector matrix V is orthogonal. That is,
+ the diagonal values of D are the eigenvalues, and
+ V*V' = I, where I is the identity matrix. The columns of V
+ represent the eigenvectors in the sense that A*V = V*D.
+
+<P>
+ If A is not symmetric, then the eigenvalue matrix D is block diagonal
+ with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
+ a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if the complex
+ eigenvalues look like
+<pre>
+
+ u + iv . . . . .
+ . u - iv . . . .
+ . . a + ib . . .
+ . . . a - ib . .
+ . . . . x .
+ . . . . . y
+</pre>
+ then D looks like
+<pre>
+
+ u v . . . .
+ -v u . . . .
+ . . a b . .
+ . . -b a . .
+ . . . . x .
+ . . . . . y
+</pre>
+ This keeps V a real matrix in both symmetric and non-symmetric
+ cases, and A*V = V*D.
+
+
+
+ <p>
+ The matrix V may be badly
+ conditioned, or even singular, so the validity of the equation
+ A = V*D*inverse(V) depends upon the condition number of V.
+
+ <p>
+ (Adapted from JAMA, a Java Matrix Library, developed by jointly
+ by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
+**/
+
+template <class Real>
+class Eigenvalue
+{
+
+
+ /** Row and column dimension (square matrix). */
+ int n;
+
+ int issymmetric; /* boolean*/
+
+ /** Arrays for internal storage of eigenvalues. */
+
+ TNT::Array1D<Real> d; /* real part */
+ TNT::Array1D<Real> e; /* img part */
+
+ /** Array for internal storage of eigenvectors. */
+ TNT::Array2D<Real> V;
+
+ /** Array for internal storage of nonsymmetric Hessenberg form.
+ @serial internal storage of nonsymmetric Hessenberg form.
+ */
+ TNT::Array2D<Real> H;
+
+
+ /** Working storage for nonsymmetric algorithm.
+ @serial working storage for nonsymmetric algorithm.
+ */
+ TNT::Array1D<Real> ort;
+
+
+ // Symmetric Householder reduction to tridiagonal form.
+
+ void tred2() {
+
+ // This is derived from the Algol procedures tred2 by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ for (int j = 0; j < n; j++) {
+ d[j] = V[n-1][j];
+ }
+
+ // Householder reduction to tridiagonal form.
+
+ for (int i = n-1; i > 0; i--) {
+
+ // Scale to avoid under/overflow.
+
+ Real scale = 0.0;
+ Real h = 0.0;
+ for (int k = 0; k < i; k++) {
+ scale = scale + fabs(d[k]);
+ }
+ if (scale == 0.0) {
+ e[i] = d[i-1];
+ for (int j = 0; j < i; j++) {
+ d[j] = V[i-1][j];
+ V[i][j] = 0.0;
+ V[j][i] = 0.0;
+ }
+ } else {
+
+ // Generate Householder vector.
+
+ for (int k = 0; k < i; k++) {
+ d[k] /= scale;
+ h += d[k] * d[k];
+ }
+ Real f = d[i-1];
+ Real g = sqrt(h);
+ if (f > 0) {
+ g = -g;
+ }
+ e[i] = scale * g;
+ h = h - f * g;
+ d[i-1] = f - g;
+ for (int j = 0; j < i; j++) {
+ e[j] = 0.0;
+ }
+
+ // Apply similarity transformation to remaining columns.
+
+ for (int j = 0; j < i; j++) {
+ f = d[j];
+ V[j][i] = f;
+ g = e[j] + V[j][j] * f;
+ for (int k = j+1; k <= i-1; k++) {
+ g += V[k][j] * d[k];
+ e[k] += V[k][j] * f;
+ }
+ e[j] = g;
+ }
+ f = 0.0;
+ for (int j = 0; j < i; j++) {
+ e[j] /= h;
+ f += e[j] * d[j];
+ }
+ Real hh = f / (h + h);
+ for (int j = 0; j < i; j++) {
+ e[j] -= hh * d[j];
+ }
+ for (int j = 0; j < i; j++) {
+ f = d[j];
+ g = e[j];
+ for (int k = j; k <= i-1; k++) {
+ V[k][j] -= (f * e[k] + g * d[k]);
+ }
+ d[j] = V[i-1][j];
+ V[i][j] = 0.0;
+ }
+ }
+ d[i] = h;
+ }
+
+ // Accumulate transformations.
+
+ for (int i = 0; i < n-1; i++) {
+ V[n-1][i] = V[i][i];
+ V[i][i] = 1.0;
+ Real h = d[i+1];
+ if (h != 0.0) {
+ for (int k = 0; k <= i; k++) {
+ d[k] = V[k][i+1] / h;
+ }
+ for (int j = 0; j <= i; j++) {
+ Real g = 0.0;
+ for (int k = 0; k <= i; k++) {
+ g += V[k][i+1] * V[k][j];
+ }
+ for (int k = 0; k <= i; k++) {
+ V[k][j] -= g * d[k];
+ }
+ }
+ }
+ for (int k = 0; k <= i; k++) {
+ V[k][i+1] = 0.0;
+ }
+ }
+ for (int j = 0; j < n; j++) {
+ d[j] = V[n-1][j];
+ V[n-1][j] = 0.0;
+ }
+ V[n-1][n-1] = 1.0;
+ e[0] = 0.0;
+ }
+
+ // Symmetric tridiagonal QL algorithm.
+
+ void tql2 () {
+
+ // This is derived from the Algol procedures tql2, by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ for (int i = 1; i < n; i++) {
+ e[i-1] = e[i];
+ }
+ e[n-1] = 0.0;
+
+ Real f = 0.0;
+ Real tst1 = 0.0;
+ Real eps = pow(2.0,-52.0);
+ for (int l = 0; l < n; l++) {
+
+ // Find small subdiagonal element
+
+ tst1 = max(tst1,fabs(d[l]) + fabs(e[l]));
+ int m = l;
+
+ // Original while-loop from Java code
+ while (m < n) {
+ if (fabs(e[m]) <= eps*tst1) {
+ break;
+ }
+ m++;
+ }
+
+
+ // If m == l, d[l] is an eigenvalue,
+ // otherwise, iterate.
+
+ if (m > l) {
+ int iter = 0;
+ do {
+ iter = iter + 1; // (Could check iteration count here.)
+
+ // Compute implicit shift
+
+ Real g = d[l];
+ Real p = (d[l+1] - g) / (2.0 * e[l]);
+ Real r = hypot(p,1.0);
+ if (p < 0) {
+ r = -r;
+ }
+ d[l] = e[l] / (p + r);
+ d[l+1] = e[l] * (p + r);
+ Real dl1 = d[l+1];
+ Real h = g - d[l];
+ for (int i = l+2; i < n; i++) {
+ d[i] -= h;
+ }
+ f = f + h;
+
+ // Implicit QL transformation.
+
+ p = d[m];
+ Real c = 1.0;
+ Real c2 = c;
+ Real c3 = c;
+ Real el1 = e[l+1];
+ Real s = 0.0;
+ Real s2 = 0.0;
+ for (int i = m-1; i >= l; i--) {
+ c3 = c2;
+ c2 = c;
+ s2 = s;
+ g = c * e[i];
+ h = c * p;
+ r = hypot(p,e[i]);
+ e[i+1] = s * r;
+ s = e[i] / r;
+ c = p / r;
+ p = c * d[i] - s * g;
+ d[i+1] = h + s * (c * g + s * d[i]);
+
+ // Accumulate transformation.
+
+ for (int k = 0; k < n; k++) {
+ h = V[k][i+1];
+ V[k][i+1] = s * V[k][i] + c * h;
+ V[k][i] = c * V[k][i] - s * h;
+ }
+ }
+ p = -s * s2 * c3 * el1 * e[l] / dl1;
+ e[l] = s * p;
+ d[l] = c * p;
+
+ // Check for convergence.
+
+ } while (fabs(e[l]) > eps*tst1);
+ }
+ d[l] = d[l] + f;
+ e[l] = 0.0;
+ }
+
+ // Sort eigenvalues and corresponding vectors.
+
+ for (int i = 0; i < n-1; i++) {
+ int k = i;
+ Real p = d[i];
+ for (int j = i+1; j < n; j++) {
+ if (d[j] < p) {
+ k = j;
+ p = d[j];
+ }
+ }
+ if (k != i) {
+ d[k] = d[i];
+ d[i] = p;
+ for (int j = 0; j < n; j++) {
+ p = V[j][i];
+ V[j][i] = V[j][k];
+ V[j][k] = p;
+ }
+ }
+ }
+ }
+
+ // Nonsymmetric reduction to Hessenberg form.
+
+ void orthes () {
+
+ // This is derived from the Algol procedures orthes and ortran,
+ // by Martin and Wilkinson, Handbook for Auto. Comp.,
+ // Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutines in EISPACK.
+
+ int low = 0;
+ int high = n-1;
+
+ for (int m = low+1; m <= high-1; m++) {
+
+ // Scale column.
+
+ Real scale = 0.0;
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list