|
|
Ublas : |
Subject: Re: [ublas] axpy_prod Weirdness
From: Vardan Akopian (vakopian_at_[hidden])
Date: 2010-08-23 18:11:00
On Mon, Aug 23, 2010 at 2:32 PM, Jane Shevtsov <jane.eco_at_[hidden]> wrote:
> 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?
>
You should not be passing the same matrix as the first and third
argument to axpy_prod, since axpy_prod uses the third argument as
output. For example with init = true, it will start by clearing your
Gprod matrix first, thus the result you're seeing.
One way to fix would be:
matrix<double> tmp(matSize, matSize);
axpy_prod(Gprod, Garr[i+j], tmp, true);
Gprod.swap(tmp);