$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/* Copyright 2008-2011 Research Foundation State University of New York   */

/* This file is part of QUB Express.                                      */

/* QUB Express is free software; you can redistribute it and/or modify    */
/* it under the terms of the GNU General Public License as published by   */
/* the Free Software Foundation, either version 3 of the License, or      */
/* (at your option) any later version.                                    */

/* QUB Express is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/* GNU General Public License for more details.                           */

/* You should have received a copy of the GNU General Public License,     */
/* named LICENSE.txt, in the QUB Express program directory.  If not, see  */
/* <http://www.gnu.org/licenses/>.                                        */

#ifdef _WIN32
  #include <windows.h>
#else
  #include <stdlib.h>
  #define BOOL int
  #define TRUE 1
  #define FALSE 0
#endif

#include <iostream>

#include "maxill_misc.h"
#include "../qubfast/ublas_matrixutil.h"

extern "C" MAXILL_API void FirstLatencies(double** A, double* pEq, double dt, int n, double** FL)
{
	// from http://www.icsu-scope.org/downloadpubs/scope34/contents.html
	//              Z = inv(I - A + col(1)*row(pEq))
	//       FL[i][j] = (Z[j][j] - Z[i][j]) / pEq[j]

	fq::matrix<double> M(n, n);
	for (int i=0; i<n; ++i)
		for (int j=0; j<n; ++j)
			M[i][j] = ((i==j)?1.0:0.0) - A[i][j] + pEq[j];
	
	fq::matrix<double> Z(n, n);
	sv_invm(M, n, Z);
	//gaussj0(M, n, NULL, 0, std::cerr);
	//for (int i=0; i<n; ++i)
	//	for (int j=0; j<n; ++j)
	//		Z[i][j] = M[i][j];

	for (int i=0; i<n; ++i)
		for (int j=0; j<n; ++j)
			FL[i][j] = dt * (Z[j][j] - Z[i][j]) / pEq[j];
}