I have a little program that implements an ecological network analysis algorithm. (It looks at the importance of indirect energy or matter flows between species.) 01 //dea.cpp -- The fundamental DEA algorithm 02 03 #include <boost/numeric/ublas/matrix.hpp> 04 #include <boost/numeric/ublas/operation.hpp> 05 #include <boost/numeric/ublas/io.hpp> //for test 06 #include "storage_adaptors.hpp" //for test 07 08 using boost::numeric::ublas::matrix; 09 10 void dea(const matrix<double>[], int, int, matrix<double>[]); 11 matrix<double> indflow(const matrix<double> & Ndea, const matrix<double> & Gdea); 12 13 int main() { 14 using boost::numeric::ublas::make_matrix_from_pointer; 15 matrix<double> Garr[4]; 16 matrix<double> Narr[2]; 17 double G1[3][3] = {0, 0, 0, 0.33333, 0, 0, 0.66667, 1, 0}; 18 double G2[3][3] = {0, 0, 0, 0.6, 0, 0, 0.4, 1, 0}; 19 double G3[3][3] = {0, 0, 0, 0.55, 0, 0, 0.45, 1, 0}; 20 double G4[3][3] = {0, 0, 0, 0.5, 0, 0, 0.5, 1, 0}; 21 Garr[0] = make_matrix_from_pointer(G1); 22 Garr[1] = make_matrix_from_pointer(G2); 23 Garr[2] = make_matrix_from_pointer(G3); 24 Garr[3] = make_matrix_from_pointer(G4); 25 dea(Garr, 2, 4, Narr); 26 std::cout << Narr[0] << std::endl; 27 matrix<double> ratio = indflow(Narr[0], Garr[0]); 28 std::cout << ratio << std::endl; 29 return 0; 30 } 31 32 //The fundamental DEA algorithm 33 //Takes array of G matrices, window size, array length and empty array for N_DEA matrices as input 34 void dea(const matrix<double> Garr[], int win, int arrLen, matrix<double> Narr[]) { 35 using namespace boost::numeric::ublas; 36 int matSize = Garr[0].size1(); 37 for (int i = 0; i < arrLen - win; i++) { 38 matrix<double> Gprod = identity_matrix<double>(matSize); 39 Narr[i] = identity_matrix<double>(matSize); 40 //the core loop 41 for (int j = 0; j < win; j++) { 42 Gprod = prod(Gprod, Garr[i+j]); 43 Narr[i] += Gprod; 44 } 45 } 46 } 47 48 //computes indirect flow fraction 49 matrix<double> indflow(const matrix<double> & Ndea, const matrix<double> & Gdea) { 50 using namespace boost::numeric::ublas; 51 int matSize = Ndea.size1(); 52 matrix<double> iflow = Ndea - Gdea - identity_matrix<double>(matSize); //Ind = N-G-I 53 matrix<double> frac = zero_matrix<double>(matSize, matSize); //inialize indirect flow fraction matrix 54 //perform division for nonzero, off-diagonal elements 55 for (int i=0; i<matSize; i++) { 56 for (int j=0; j<matSize; j++) { 57 if (i != j && Ndea(i,j) != 0) //avoid division by zero 58 frac(i, j) = iflow(i, j) / Ndea(i, j); 59 } 60 } 61 return frac; 62 } The correct output is below. [3,3]((1,0,0),(0.33333,1,0),(1.26667,1,1)) [3,3]((0,0,0),(0,0,0),(0.473683,0,0)) The code given above works, but to speed it up, I tried replacing the first line of the core loop in dea() with axpy_prod(Gprod, Garr[i+j], Gprod, true);. This gives the following wrong output: [3,3]((1,0,0),(0,1,0),(0,0,1)) [3,3]((0,0,0),(0,0,0),(0,0,0)) And if I use axpy_prod(Gprod, Garr[i+j], Gprod, false);, I get: [3,3]((3,0,0),(1.26666,3,0),(2.33334,3,3)) [3,3]((0,0,0),(0.736843,0,0),(0.714285,0.666667,0)) What gives? When should I avoid using axpy_prod? -- ------------- Jane Shevtsov Ecology Ph.D. candidate, University of Georgia co-founder, <www.worldbeyondborders.org> Check out my blog, <http://perceivingwholes.blogspot.com>Perceiving Wholes "The whole person must have both the humility to nurture the Earth and the pride to go to Mars." --Wyn Wachhorst, The Dream of Spaceflight