[Boost-users] axpy_prod Weirdness

Subject: [Boost-users] axpy_prod Weirdness
From: Jane Shevtsov (jane.eco_at_[hidden])
Date: 2010-08-22 15:13:10

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