directsolve.cpp

斯坦福Energy211/CME211课《c++编程——地球科学科学家和工程师》的课件

#include <stdexcept>	// for exception handling
#include "mex.h"
#include "linsys.h"

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
	// Get size of matrix passed as first argument
	int m = mxGetM( prhs[0] );
	int n = mxGetN( prhs[0] );
	
	// Get pointers to matrix and rhs data
	double *pA = mxGetPr( prhs[0] );
	double *pb = mxGetPr( prhs[1] );
	
	// Matrices in MATLAB are stored in COLUMN order
	Matrix A( m, n );
	for ( int i = 0; i < m; i++ )
		for ( int j = 0; j < n; j++ )
			A[i][j] = pA[j * m + i];
	ColVector b( m, pb );
	ColVector x( n );
	
	try
	{
		LinearSystem S( A );
		x = S.DirectSolve( b );
	
		// Create MATLAB vector for output argument
		plhs[0] = mxCreateDoubleMatrix( n, 1, mxREAL );
		double *px = mxGetPr( plhs[0] );
		for ( int i = 0; i < n; i++ )
			px[i] = x[i];
	}
	catch( std::runtime_error e )
	{
		// Pass error message up to MATLAB
		mexErrMsgTxt( e.what() );
	}
}