[MOAB-dev] commit/MOAB: iulian07: include EigenDecomp into Matrix3 and move Push F90 example into imesh
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Nov 27 10:30:24 CST 2013
1 new commit in MOAB:
https://bitbucket.org/fathomteam/moab/commits/095efb6f6939/
Changeset: 095efb6f6939
Branch: master
User: iulian07
Date: 2013-11-27 17:24:48
Summary: include EigenDecomp into Matrix3 and move Push F90 example into imesh
so remove EigenDecomp.hpp from the repostory, and
move parallel fortran90 example into itaps/imesh
The advantage is that the autotools will figure out the correct
compiling and linking
the parallel test can be done after compiling (or make check) as before,
mpiexec -np 2 Push...
For the time being, the file is still in example folder
The compiling instructions work on mpi+gcc, but intel compiler on fusion was
still unhappy.
Affected #: 6 files
diff --git a/itaps/imesh/Makefile.am b/itaps/imesh/Makefile.am
index 647b8c1..d8a81eb 100644
--- a/itaps/imesh/Makefile.am
+++ b/itaps/imesh/Makefile.am
@@ -38,6 +38,9 @@ if PARALLEL
# check_PROGRAMS += partest iMeshP_unit_tests moabtest
check_PROGRAMS += partest MOAB_iMeshP_unit_tests
+if ENABLE_FORTRAN
+ check_PROGRAMS += PushParMeshIntoMoabF90
+endif
# moabtest_SOURCES = moabtest.F
# check_PROGRAMS += ftest
@@ -69,6 +72,9 @@ ScdMeshF77_DEPENDENCIES = $(TESTDEPS)
FindAdjacencyF90_SOURCES = FindAdjacencyF90.F90
FindAdjacencyF90_DEPENDENCIES = $(TESTDEPS)
+PushParMeshIntoMoabF90_SOURCES = PushParMeshIntoMoabF90.F90
+PushParMeshIntoMoabF90_DEPENDENCIES = $(TESTDEPS)
+
lib_LTLIBRARIES = libiMesh.la
libiMesh_la_includedir = $(includedir)
diff --git a/itaps/imesh/PushParMeshIntoMoabF90.F90 b/itaps/imesh/PushParMeshIntoMoabF90.F90
new file mode 100644
index 0000000..46d577a
--- /dev/null
+++ b/itaps/imesh/PushParMeshIntoMoabF90.F90
@@ -0,0 +1,283 @@
+! PushParMeshIntoMoabF90: push parallel mesh into moab, F90 version
+!
+! This program shows how to push a mesh into MOAB in parallel from Fortran90, with sufficient
+! information to resolve boundary sharing and exchange a layer of ghost information.
+! To successfully link this example, you need to specify FCFLAGS that include:
+! a) -DUSE_MPI, and
+! b) flags required to link Fortran90 MPI programs with the C++ compiler; these flags
+! can often be found on your system by inspecting the output of 'mpif90 -show'
+! For example, using gcc, the link line looks like:
+! make MOAB_DIR=<moab install dir> FCFLAGS="-DUSE_MPI -I/usr/lib/openmpi/include -pthread -I/usr/lib/openmpi/lib -L/usr/lib/openmpi/lib -lmpi_f90 -lmpi_f77 -lmpi -lopen-rte -lopen-pal -ldl -Wl,--export-dynamic -lnsl -lutil -lm -ldl" PushParMeshIntoMoabF90
+!
+! Usage: PushParMeshIntoMoab
+
+program PushParMeshIntoMoab
+
+ use ISO_C_BINDING
+ implicit none
+
+#include "mpif.h"
+
+#ifdef USE_MPI
+# include "iMeshP_f.h"
+#else
+# include "iMesh_f.h"
+#endif
+
+ ! declarations
+ ! imesh is the instance handle
+ iMesh_Instance imesh
+ ! NUMV, NUME, NVPERE are the hardwired here; these are for the whole mesh,
+ ! local mesh determined later
+ integer NUMV, NUME, NVPERE
+ parameter (NUMV = 8) ! # vertices in whole mesh
+ parameter (NUME = 6) ! # elements in whole mesh
+ parameter (NVPERE = 4) ! # vertices per element
+ ! ents, verts will be arrays storing vertex/entity handles
+ iBase_EntityHandle, pointer :: ents, verts
+ iBase_EntitySetHandle root_set
+ TYPE(C_PTR) :: vertsPtr, entsPtr
+ ! storage for vertex positions, element connectivity indices, global vertex ids
+ real*8 coords(0:3*NUMV-1)
+ integer iconn(0:4*NUME-1), gids(0:NUMV-1)
+ !
+ ! local variables
+ integer lgids(0:NUMV-1), lconn(0:4*NUME-1)
+ real*8 lcoords(0:3*NUMV-1)
+ integer lnumv, lvids(0:NUMV-1), gvids(0:NUMV-1)
+ integer lvpe, ltp ! lvpe = # vertices per entity, ltp = element type
+ integer ic, ie, iv, istart, iend, ierr, indv, lnume, rank, sz
+
+#ifdef USE_MPI
+ ! local variables for parallel runs
+ iMeshP_PartitionHandle imeshp
+ ! integer MPI_COMM_WORLD
+#endif
+
+ ! vertex positions, latlon coords, (lat, lon, lev), fortran ordering
+ ! (first index varying fastest)
+ data coords / &
+ 0.0, -45.0, 0.0, 90.0, -45.0, 0.0, 180.0, -45.0, 0.0, 270.0, -45.0, 0.0, &
+ 0.0, 45.0, 0.0, 90.0, 45.0, 0.0, 180.0, 45.0, 0.0, 270.0, 45.0, 0.0 /
+
+ ! quad index numbering, each quad ccw, sides then bottom then top
+ data iconn / &
+ 0, 1, 5, 4, &
+ 1, 2, 6, 5, &
+ 2, 3, 7, 6, &
+ 3, 0, 4, 7, &
+ 0, 3, 2, 1, &
+ 4, 5, 6, 7 /
+
+ data lvpe /4/ ! quads in this example
+ data ltp / iMesh_QUADRILATERAL / ! from iBase_f.h
+
+ ! initialize global vertex ids
+ do iv = 0, NUMV-1
+ lgids(iv) = iv+1
+ end do
+
+#ifdef USE_MPI
+ ! init the parallel partition
+ call MPI_INIT(ierr)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD, sz, ierr)
+ call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
+ ! compute starting/ending element numbers
+ lnume = NUME / sz
+ istart = rank * lnume
+ iend = istart + lnume - 1
+ if (rank .eq. sz-1) then
+ iend = NUME-1
+ lnume = iend - istart + 1
+ endif
+#else
+ ! set the starting/ending element numbers
+ istart = 0
+ iend = NUME-1
+ lnume = NUME
+#endif
+
+ ! for my elements, figure out which vertices I use and accumulate local indices and coords
+ ! lvids stores the local 0-based index for each vertex; -1 means vertex i isn't used locally
+ ! also build up connectivity indices for local elements, in lconn
+ do iv = 0, NUMV-1
+ lvids(iv) = -1
+ end do
+ lnumv = -1
+ do ie = istart, iend
+ do iv = 0, lvpe-1
+ indv = iconn(lvpe*ie + iv)
+ if (lvids(indv) .eq. -1) then
+ lnumv = lnumv + 1 ! increment local # verts
+ do ic = 0, 2 ! cache local coords
+ lcoords(3*lnumv+ic) = coords(3*indv+ic)
+ end do
+ lvids(indv) = lnumv
+ gvids(lnumv) = 1+indv
+ end if
+ lconn(lvpe*(ie-istart)+iv) = lvids(indv)
+ end do ! do iv
+ end do ! do ie
+
+ lnumv = lnumv + 1
+
+ ! now create the mesh; this also initializes parallel sharing and ghost exchange
+ imesh = 0
+ imeshp = 0
+ call create_mesh(imesh, imeshp, MPI_COMM_WORLD, lnumv, lnume, gvids, lvpe, ltp, lcoords, lconn, &
+ vertsPtr, entsPtr, ierr)
+ call c_f_pointer(vertsPtr, verts, [lnumv])
+ call c_f_pointer(entsPtr, ents, [lnume])
+
+ ! get/report number of vertices, elements
+ call iMesh_getRootSet(%VAL(imesh), root_set, ierr)
+ iv = 0
+ ie = 0
+#ifdef USE_MPI
+ call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr)
+ call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_FACE), ie, ierr)
+ if (rank .eq. 0) then
+ write(0,*) "Number of vertices = ", iv
+ write(0,*) "Number of entities = ", ie
+ endif
+#else
+ call iMesh_getNumOfTypeAll(%VAL(imesh), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr)
+ call iMesh_getNumOfTypeAll(%VAL(imesh), %VAL(root_set), %VAL(iBase_FACE), ie, ierr)
+ write(0,*) "Number of vertices = ", iv
+ write(0,*) "Number of entities = ", ie
+#endif
+
+ ! from here, can use verts and ents as (1-based) arrays of entity handles for input to other iMesh functions
+
+ call MPI_FINALIZE(ierr)
+ stop
+end program PushParMeshIntoMoab
+
+subroutine create_mesh( &
+ ! interfaces
+ imesh, imeshp, &
+ ! input
+ comm, numv, nume, vgids, nvpe, tp, posn, iconn, &
+ ! output
+ vertsPtr, entsPtr, ierr)
+ !
+ ! create a mesh with numv vertices and nume elements, with elements of type tp
+ ! vertices have positions in posn (3 coordinates each, interleaved xyzxyz...), indexed from 0
+ ! elements have nvpe vertices per entity, with connectivity indices stored in iconn, referencing
+ ! vertices using 0-based indices; vertex and entity handles are output in arrays passed in
+ !
+ ! if imesh/imeshp are 0, imesh/imeshp are initialized in this subroutine
+ !
+
+ use ISO_C_BINDING
+ implicit none
+
+#ifdef USE_MPI
+# include "iMeshP_f.h"
+# include "mpif.h"
+#else
+# include "iMesh_f.h"
+#endif
+
+ ! subroutine arguments
+ iMesh_Instance imesh
+ TYPE(C_PTR) :: vertsPtr, entsPtr
+ integer numv, nume, nvpe, vgids(0:*), iconn(0:*), ierr, tp
+ real*8 posn(0:*)
+#ifdef USE_MPI
+ iMeshP_PartitionHandle imeshp
+ integer comm
+#endif
+
+ ! local variables
+ integer comm_sz, comm_rank, numa, numo, iv, ie
+ TYPE(C_PTR) :: statsPtr
+ integer, allocatable, target :: stats(:)
+ iBase_TagHandle tagh
+ integer i
+ iBase_EntityHandle, pointer :: verts(:), ents(:)
+ iBase_EntityHandle, allocatable :: conn(:)
+ iBase_EntitySetHandle root_set
+#ifdef USE_MPI
+ IBASE_HANDLE_T mpi_comm_c
+ TYPE(C_PTR) :: partsPtr
+ iMeshP_PartHandle, pointer :: parts(:)
+ iMeshP_PartHandle part
+ integer partsa, partso
+#endif
+
+ ! create the Mesh instance
+ if (imesh .eq. 0) then
+ call iMesh_newMesh("MOAB", imesh, ierr)
+ end if
+
+#ifdef USE_MPI
+ if (imeshp .eq. 0) then
+ call iMeshP_getCommunicator(%VAL(imesh), MPI_COMM_WORLD, mpi_comm_c, ierr)
+ call iMeshP_createPartitionAll(%VAL(imesh), %VAL(mpi_comm_c), imeshp, ierr)
+ call iMeshP_createPart(%VAL(imesh), %VAL(imeshp), part, ierr)
+ else
+ partsa = 0
+ call iMeshP_getLocalParts(%VAL(imesh), %VAL(imeshp), partsPtr, partsa, partso, ierr)
+ call c_f_pointer(partsPtr, parts, [partso])
+ part = parts(1)
+ end if
+ call MPI_COMM_RANK(comm, comm_rank, ierr)
+ call MPI_COMM_SIZE(comm, comm_sz, ierr)
+#endif
+
+ ! create the vertices, all in one call
+ numa = 0
+ call iMesh_createVtxArr(%VAL(imesh), %VAL(numv), %VAL(iBase_INTERLEAVED), posn, %VAL(3*numv), &
+ vertsPtr, numa, numo, ierr)
+
+ ! fill in the connectivity array, based on indexing from iconn
+ allocate (conn(0:nvpe*nume-1))
+ call c_f_pointer(vertsPtr, verts, [numv])
+ do i = 0, nvpe*nume-1
+ conn(i) = verts(1+iconn(i))
+ end do
+ ! create the elements
+ numa = 0
+ allocate(stats(0:nume-1))
+ statsPtr = C_LOC(stats(0))
+ call iMesh_createEntArr(%VAL(imesh), %VAL(tp), conn, %VAL(nvpe*nume), &
+ entsPtr, numa, numo, statsPtr, numa, numo, ierr)
+ deallocate(stats)
+ deallocate(conn)
+
+#ifdef USE_MPI
+ ! take care of parallel stuff
+
+ ! add entities to part, using iMesh
+ call c_f_pointer(entsPtr, ents, [numo])
+ call iMesh_addEntArrToSet(%VAL(imesh), ents, %VAL(numo), %VAL(part), ierr)
+ ! set global ids on vertices, needed for sharing between procs
+ call iMesh_getTagHandle(%VAL(imesh), "GLOBAL_ID", tagh, ierr, %VAL(9))
+ if (iBase_SUCCESS .ne. ierr) then
+ ! didn't get handle, need to create the tag
+ call iMesh_createTag(%VAL(imesh), "GLOBAL_ID", %VAL(iBase_INTEGER), tagh, ierr, %VAL(9))
+ end if
+ call iMesh_setIntArrData(%VAL(imesh), verts, %VAL(numv), %VAL(tagh), vgids, %VAL(numv), ierr)
+
+ ! now resolve shared verts and exchange ghost cells
+ call iMeshP_syncMeshAll(%VAL(imesh), %VAL(imeshp), ierr)
+ call iMesh_getRootSet(%VAL(imesh), root_set, ierr)
+ call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_VERTEX), iv, ierr)
+ call iMeshP_getNumOfTypeAll(%VAL(imesh), %VAL(imeshp), %VAL(root_set), %VAL(iBase_FACE), ie, ierr)
+ if (comm_rank .eq. 0) then
+ write(0,*) "After syncMeshAll:"
+ write(0,*) " Number of vertices = ", iv
+ write(0,*) " Number of entities = ", ie
+ endif
+
+ call iMeshP_createGhostEntsAll(%VAL(imesh), %VAL(imeshp), %VAL(2), %VAL(1), %VAL(1), %VAL(0), ierr)
+ if (comm_rank .eq. 0) then
+ write(0,*) "After createGhostEntsAll:"
+ write(0,*) " Number of vertices = ", iv
+ write(0,*) " Number of entities = ", ie
+ endif
+#endif
+
+ return
+end subroutine create_mesh
diff --git a/src/Makefile.am b/src/Makefile.am
index 64b8a0f..0eb07c5 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -157,7 +157,6 @@ nobase_libMOAB_la_include_HEADERS = \
moab/GeomTopoTool.hpp \
moab/HigherOrderFactory.hpp \
moab/HomXform.hpp \
- moab/EigenDecomp.hpp \
moab/EntityType.hpp \
moab/EntityHandle.hpp \
moab/FBEngine.hpp \
diff --git a/src/OrientedBox.cpp b/src/OrientedBox.cpp
index caab7d1..0932952 100644
--- a/src/OrientedBox.cpp
+++ b/src/OrientedBox.cpp
@@ -274,7 +274,7 @@ ErrorCode OrientedBox::compute_from_vertices( OrientedBox& result,
// Get axes (Eigenvectors) from covariance matrix
double lambda[3];
- EigenDecomp( a, lambda, result.axis );
+ moab::Matrix::EigenDecomp( a, lambda, result.axis );
// Calculate center and extents of box given orientation defined by axes
return box_from_axes( result, instance, vertices );
@@ -369,7 +369,7 @@ ErrorCode OrientedBox::compute_from_covariance_data(
// get axes (Eigenvectors) from covariance matrix
double lamda[3];
- EigenDecomp( data.matrix, lamda, result.axis );
+ moab::Matrix::EigenDecomp( data.matrix, lamda, result.axis );
// We now have only the axes. Calculate proper center
// and extents for enclosed points.
diff --git a/src/moab/EigenDecomp.hpp b/src/moab/EigenDecomp.hpp
deleted file mode 100644
index 12110eb..0000000
--- a/src/moab/EigenDecomp.hpp
+++ /dev/null
@@ -1,165 +0,0 @@
-/*
- * MOAB, a Mesh-Oriented datABase, is a software component for creating,
- * storing and accessing finite element mesh data.
- *
- * Copyright 2004 Sandia Corporation. Under the terms of Contract
- * DE-AC04-94AL85000 with Sandia Coroporation, the U.S. Government
- * retains certain rights in this software.
- *
- * This library is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public
- * License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- *
- */
-
-#ifndef MB_EIGENDECOMP_HPP
-#define MB_EIGENDECOMP_HPP
-
-#include <math.h>
-#include <iostream>
-
-namespace moab {
-//TODO: is this copied code actually able to be licensed
-//TODO: Actually should have license text here, not link
-// Copied from "Visualization Toolkit"
-// Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
-// All rights reserved.
-// See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
-//
-// Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
-// real symmetric matrix. Square nxn matrix a; size of matrix in n;
-// output eigenvalues in w; and output eigenvectors in v. Resulting
-// eigenvalues/vectors are sorted in decreasing order; eigenvectors are
-// normalized.
-// TODO: Remove this
-#define VTK_ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
- a[k][l]=h+s*(g-h*tau)
-
-//TODO: Refactor this method into subroutines
-//use a namespace { } with no name to
-//contain subroutines so that the compiler
-//automatically inlines them.
-
-template< typename Matrix, typename Vector>
-ErrorCode EigenDecomp( const Matrix & _a,
- double w[3],
- Vector v[3] ) {
- const int MAX_ROTATIONS = 20;
- const double one_ninth = 1./9;
- int i, j, k, iq, ip, numPos;
- double tresh, theta, tau, t, sm, s, h, g, c, tmp;
- double b[3], z[3];
- Matrix a( _a);
-
- // initialize
- for (ip=0; ip<3; ip++) {
- for (iq=0; iq<3; iq++){
- v[ip][iq] = 0.0;
- }
- v[ip][ip] = 1.0;
- }
- for (ip=0; ip<3; ip++) {
- b[ip] = w[ip] = a[ip][ip];
- z[ip] = 0.0;
- }
-
- // begin rotation sequence
- for (i=0; i<MAX_ROTATIONS; i++){
- sm = 0.0;
- for (ip=0; ip<2; ip++){
- for (iq=ip+1; iq<3; iq++){ sm += fabs(a[ip][iq]); }
- }
-
- if ( sm == 0.0 ){ break; }
- // first 3 sweeps
- tresh = (i < 3)? 0.2*sm*one_ninth : 0.0;
- for (ip=0; ip<2; ip++) {
- for (iq=ip+1; iq<3; iq++) {
- g = 100.0*fabs(a[ip][iq]);
-
- // after 4 sweeps
- if ( i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
- && (fabs(w[iq])+g) == fabs(w[iq])) {
- a[ip][iq] = 0.0;
- }
- else if ( fabs(a[ip][iq]) > tresh) {
- h = w[iq] - w[ip];
- if ( (fabs(h)+g) == fabs(h)){ t = (a[ip][iq]) / h; }
- else {
- theta = 0.5*h / (a[ip][iq]);
- t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
- if (theta < 0.0) { t = -t;}
- }
- c = 1.0 / sqrt(1+t*t);
- s = t*c;
- tau = s/(1.0+c);
- h = t*a[ip][iq];
- z[ip] -= h;
- z[iq] += h;
- w[ip] -= h;
- w[iq] += h;
- a[ip][iq]=0.0;
- // ip already shifted left by 1 unit
- for (j = 0;j <= ip-1;j++) { VTK_ROTATE(a,j,ip,j,iq); }
- // ip and iq already shifted left by 1 unit
- for (j = ip+1;j <= iq-1;j++) { VTK_ROTATE(a,ip,j,j,iq); }
- // iq already shifted left by 1 unit
- for (j=iq+1; j<3; j++) { VTK_ROTATE(a,ip,j,iq,j); }
- for (j=0; j<3; j++) { VTK_ROTATE(v,j,ip,j,iq); }
- }
- }
- }
-
- for (ip=0; ip<3; ip++) {
- b[ip] += z[ip];
- w[ip] = b[ip];
- z[ip] = 0.0;
- }
- }
-
- //// this is NEVER called
- if ( i >= MAX_ROTATIONS ) {
- std::cerr << "Matrix3D: Error extracting eigenfunctions" << std::endl;
- return MB_FAILURE;
- }
-
- // sort eigenfunctions these changes do not affect accuracy
- for (j=0; j<2; j++){ // boundary incorrect
- k = j;
- tmp = w[k];
- for (i=j+1; i<3; i++){ // boundary incorrect, shifted already
- if (w[i] >= tmp){ // why exchage if same?
- k = i;
- tmp = w[k];
- }
- }
- if (k != j){
- w[k] = w[j];
- w[j] = tmp;
- for (i=0; i<3; i++){
- tmp = v[i][j];
- v[i][j] = v[i][k];
- v[i][k] = tmp;
- }
- }
- }
- // insure eigenvector consistency (i.e., Jacobi can compute vectors that
- // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
- // reek havoc in hyperstreamline/other stuff. We will select the most
- // positive eigenvector.
- int ceil_half_n = (3 >> 1) + (3 & 1);
- for (j=0; j<3; j++) {
- for (numPos=0, i=0; i<3; i++) {
- if ( v[i][j] >= 0.0 ) { numPos++; }
- }
-// if ( numPos < ceil(double(n)/double(2.0)) )
- if ( numPos < ceil_half_n) {
- for(i=0; i<3; i++) { v[i][j] *= -1.0; }
- }
- }
- return MB_SUCCESS;
-}
-
-} // namespace moab
-#endif //MB_EIGENDECOMP_HPP
diff --git a/src/moab/Matrix3.hpp b/src/moab/Matrix3.hpp
index 08d38da..0272159 100644
--- a/src/moab/Matrix3.hpp
+++ b/src/moab/Matrix3.hpp
@@ -24,7 +24,8 @@
#define MOAB_MATRIX3_HPP
#include "moab/Types.hpp"
-#include "moab/EigenDecomp.hpp"
+//#include "moab/EigenDecomp.hpp"
+#include <iostream>
#include "moab/CartVect.hpp"
#include <iosfwd>
@@ -121,6 +122,142 @@ namespace Matrix{
res[ 2] = v[0] * m(2,0) + v[1] * m(2,1) + v[2] * m(2,2);
return res;
}
+
+ // moved from EigenDecomp.hpp
+
+ // Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
+ // real symmetric matrix. Square nxn matrix a; size of matrix in n;
+ // output eigenvalues in w; and output eigenvectors in v. Resulting
+ // eigenvalues/vectors are sorted in decreasing order; eigenvectors are
+ // normalized.
+ // TODO: Remove this
+ #define VTK_ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
+ a[k][l]=h+s*(g-h*tau)
+
+ //TODO: Refactor this method into subroutines
+ //use a namespace { } with no name to
+ //contain subroutines so that the compiler
+ //automatically inlines them.
+
+ template< typename Matrix, typename Vector>
+ ErrorCode EigenDecomp( const Matrix & _a,
+ double w[3],
+ Vector v[3] ) {
+ const int MAX_ROTATIONS = 20;
+ const double one_ninth = 1./9;
+ int i, j, k, iq, ip, numPos;
+ double tresh, theta, tau, t, sm, s, h, g, c, tmp;
+ double b[3], z[3];
+ Matrix a( _a);
+
+ // initialize
+ for (ip=0; ip<3; ip++) {
+ for (iq=0; iq<3; iq++){
+ v[ip][iq] = 0.0;
+ }
+ v[ip][ip] = 1.0;
+ }
+ for (ip=0; ip<3; ip++) {
+ b[ip] = w[ip] = a[ip][ip];
+ z[ip] = 0.0;
+ }
+
+ // begin rotation sequence
+ for (i=0; i<MAX_ROTATIONS; i++){
+ sm = 0.0;
+ for (ip=0; ip<2; ip++){
+ for (iq=ip+1; iq<3; iq++){ sm += fabs(a[ip][iq]); }
+ }
+
+ if ( sm == 0.0 ){ break; }
+ // first 3 sweeps
+ tresh = (i < 3)? 0.2*sm*one_ninth : 0.0;
+ for (ip=0; ip<2; ip++) {
+ for (iq=ip+1; iq<3; iq++) {
+ g = 100.0*fabs(a[ip][iq]);
+
+ // after 4 sweeps
+ if ( i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
+ && (fabs(w[iq])+g) == fabs(w[iq])) {
+ a[ip][iq] = 0.0;
+ }
+ else if ( fabs(a[ip][iq]) > tresh) {
+ h = w[iq] - w[ip];
+ if ( (fabs(h)+g) == fabs(h)){ t = (a[ip][iq]) / h; }
+ else {
+ theta = 0.5*h / (a[ip][iq]);
+ t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
+ if (theta < 0.0) { t = -t;}
+ }
+ c = 1.0 / sqrt(1+t*t);
+ s = t*c;
+ tau = s/(1.0+c);
+ h = t*a[ip][iq];
+ z[ip] -= h;
+ z[iq] += h;
+ w[ip] -= h;
+ w[iq] += h;
+ a[ip][iq]=0.0;
+ // ip already shifted left by 1 unit
+ for (j = 0;j <= ip-1;j++) { VTK_ROTATE(a,j,ip,j,iq); }
+ // ip and iq already shifted left by 1 unit
+ for (j = ip+1;j <= iq-1;j++) { VTK_ROTATE(a,ip,j,j,iq); }
+ // iq already shifted left by 1 unit
+ for (j=iq+1; j<3; j++) { VTK_ROTATE(a,ip,j,iq,j); }
+ for (j=0; j<3; j++) { VTK_ROTATE(v,j,ip,j,iq); }
+ }
+ }
+ }
+
+ for (ip=0; ip<3; ip++) {
+ b[ip] += z[ip];
+ w[ip] = b[ip];
+ z[ip] = 0.0;
+ }
+ }
+
+ //// this is NEVER called
+ if ( i >= MAX_ROTATIONS ) {
+ std::cerr << "Matrix3D: Error extracting eigenfunctions" << std::endl;
+ return MB_FAILURE;
+ }
+
+ // sort eigenfunctions these changes do not affect accuracy
+ for (j=0; j<2; j++){ // boundary incorrect
+ k = j;
+ tmp = w[k];
+ for (i=j+1; i<3; i++){ // boundary incorrect, shifted already
+ if (w[i] >= tmp){ // why exchage if same?
+ k = i;
+ tmp = w[k];
+ }
+ }
+ if (k != j){
+ w[k] = w[j];
+ w[j] = tmp;
+ for (i=0; i<3; i++){
+ tmp = v[i][j];
+ v[i][j] = v[i][k];
+ v[i][k] = tmp;
+ }
+ }
+ }
+ // insure eigenvector consistency (i.e., Jacobi can compute vectors that
+ // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
+ // reek havoc in hyperstreamline/other stuff. We will select the most
+ // positive eigenvector.
+ int ceil_half_n = (3 >> 1) + (3 & 1);
+ for (j=0; j<3; j++) {
+ for (numPos=0, i=0; i<3; i++) {
+ if ( v[i][j] >= 0.0 ) { numPos++; }
+ }
+ // if ( numPos < ceil(double(n)/double(2.0)) )
+ if ( numPos < ceil_half_n) {
+ for(i=0; i<3; i++) { v[i][j] *= -1.0; }
+ }
+ }
+ return MB_SUCCESS;
+ }
} //namespace Matrix
class Matrix3 {
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list