Skip to content

Instantly share code, notes, and snippets.

@jashwanth9
Created August 7, 2012 18:11
Show Gist options
  • Save jashwanth9/3287945 to your computer and use it in GitHub Desktop.
Save jashwanth9/3287945 to your computer and use it in GitHub Desktop.
namespace dspsv_func{
const int PRINT_DEBUG_STUFF = 0;
function dspsv_exec(var A,var b,int uplo)
{
var pla = loadlib("linalg_group");
var lapack = loadlib('liblapack.so');
if (lapack == null || !lapack)
die("Cannot find liblapack. Search for the correct paths");
int m,n,nrhs,ldb,info;
m=A.rows;
n=A.cols;
nrhs=b.cols;
int ipiv[n];
int j;
for(j=0;j<n;++j)
ipiv[j]=0;
ldb=max(1,n);
say(b);
var AP;
AP=create_ap(A,uplo);
say(AP);
var dspsv = dlfunc(lapack, "dspsv_", "vpppppppp");
if(dspsv == null || !dspsv)
die("Not DSPSV");
dspsv(uplo,n,nrhs,AP,ipiv,b,ldb,info);
if(info==0)
return AP;
}
function create_ap(var A,int uplo)
{
int m,n;
m=A.rows;
n=A.cols;
var pla = loadlib("linalg_group");
int ap_size;
ap_size=((n*(n+1))/2);
var ap=new 'NumMatrix2D';
ap.resize(1,ap_size);
int i,j;
if(uplo == 'U')
{
for(i=0;i<m;++i)
for(j=0;j<n;++j)
ap[i + (j-1)*j/2] = A[i,j];
}
else if(uplo == 'L')
{
for(i=0;i<m;++i)
for(j=0;j<n;++j)
ap[i + (j-1)*((2*n)-j)/2] = A[i,j];
}
else
{
die("invalid input");
}
return ap;
}
function max(var a,var b)
{
return a>b?a:b;
}
function debug(var matrix, string msg, var args [slurpy])
{
if (PRINT_DEBUG_STUFF) {
say(sprintf(msg, args));
say(matrix);
}
}
}
function main[main](var args)
{
var pla = loadlib("linalg_group");
var a = new 'NumMatrix2D';
var b=new 'NumMatrix2D';
var ap=new 'NumMatrix2D';
int uplo;
uplo=ord('U',0);
a.initialize_from_args(3, 3,
1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0);
b.initialize_from_args(3, 1,
1.0,
2.0,
3.0);
using dspsv_func.dspsv_exec;
ap=dspsv_exec(a,b);
say("successful");
say("AP=");
say(ap);
say("B=");
say(b);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment