mathematicaをLAPACK代わりに使う

一般化固有値問題をLAPACKを使って解いているのですが、行列が非対称になったときに
たまに結果が収束しなかったりして困ってました。僕のLAPACKの使い方が駄目なだけかもしれませんが、
倍精度doubleでは駄目かもしれない。ちなみに、単精度(float)は全然話になりません。
４倍精度のlong doubleなるものがあるのですが、今度はLAPACKが対応していない。
４倍精度に対応したMLAPACKや、別の固有値演算パック(Arnoldi法で解くやつ)のARPACKとか
あるようですが、FORTRANで書かれていて、僕の環境ではVCからのコンパイルやリンクの仕方が
意味不明すぎる。

ほとほと困っていたのですが、そういえば職場ではmathematicaが使えて、これで固有値問題を
解かせるとちゃんと結果が出る。そこで、これをLAPACK代わりに使えないかと調べたところ、
mathlinkなるもので簡単に使えるようです。
以下に、mathlinkを使った一般化固有値問題を解くコードを置きます。
VCではmathematicaのライブラリとインクルードファイルのパスを通しておいてください。
使い方はinit_kernelでmathematicaのカーネルを起動して、solve_generalized_eigen_problemで
カーネルに固有値問題を解くようにメッセージを送って,戻ってきた結果を引数にコピーしているだけです。
メッセージ内容を書き換えると、ほかの数値計算にも使えると思います。

//
// mathematicaを使う
//
#include <float.h>
// mathlink読み出し
#include "mathlink.h"

// グローバル変数
MLINK link;
MLENV env;
char *str_S, *str_H;
char *str_in;

// NxN行列データをmathematica形式のテクストデータに変換
void conv_mat2str(char *str, double *A, int N)
{
	int i, j, n, s;
	n = 0;
	s = sprintf(&str[n], "{");
	n += s;
	for(i = 0;i < N;i++){
		s = sprintf(&str[n], "{");
		n += s;
		for(j = 0;j < N;j++){			
			if(j != (N-1))
				s = sprintf(&str[n], "%.18f, ", A[i * N + j]);
			else
				s = sprintf(&str[n], "%.18f}",  A[i * N + j]);
			n += s;
		};//j
		if(i != (N-1))
			s = sprintf(&str[n], ",");
		else
			s = sprintf(&str[n], "}\0");
		n+= s;
	};//i
};

// カーネル起動
void init_kernel()
{
	char *in[] = {
		"-linklaunch",
			"-linkname",
			"\"C:\\Program Files\\Wolfram Research\\Mathematica\\7.0\\MathKernel\" -mathlink",
			""
	};
	int errno;
	env = MLInitialize((char *)0);
	// math kernelを開く
	link = MLOpenArgv(env, in, in + 3, &errno);
	// テンポラリテキストデータワークを適当に確保
	str_S  = new char[30000];
	str_H  = new char[30000];
	str_in = new char[30000];
};

// カーネル閉じる
void close_kernel()
{
	// mathkernelを閉じる
	MLClose(link);
	MLDeinitialize(env);
	//ワーク解放
	delete str_S;
	delete str_H;
	delete str_in;
}

// 一般化固有値問題Hx=λSxを解く(実数値のみ取り出す)
void solve_generalized_eigen_problem(double *H, double *S, int N, double *E)
{
	// 行列データをmathematicaの入力テキストデータに変換
	conv_mat2str(str_S, S, N);
	conv_mat2str(str_H, H, N);
	// methematicaに入力する式作成(※注意：Re[]で無理やり実数化してます)
	// [1] A = Re[Eigensystem[{H, S}]];
	// [2] A[[1]];固有値
	// [3] A[[2]];固有ベクトル
	sprintf(str_in, "A=Re[Eigensystem[{%s, %s}]];A[[1]]\0", str_H, str_S);
    /* math kernelに計算式を送信 */
	MLPutFunction(link, "EvaluatePacket", 1);
    MLPutFunction(link, "ToExpression", 1);
    MLPutString(link, (const char *)str_in);
    MLEndPacket(link);
    /* send the packet */
	MLFlush(link);
	// 結果を受け取る
	int p;
	while ((p = MLNextPacket(link)) && p != RETURNPKT){
		MLNewPacket(link);
	};
	// 結果を変数にコピー
	double *ret;
	long n;
	double *data;
	long *dims;
	char **heads;
	long d;
	// 固有値データを受け取る
	MLGetRealList(link, &ret, &n);
    /* math kernelに計算式を送信 */
	MLPutFunction(link, "EvaluatePacket", 1);
    MLPutFunction(link, "ToExpression", 1);
    MLPutString(link, (const char *)"A[[2]]");
    MLEndPacket(link);
    /* send the packet */
	MLFlush(link);
	// 結果を受け取る
	while ((p = MLNextPacket(link)) && p != RETURNPKT){
		MLNewPacket(link);
	};
	// 固有ベクトルデータ
	MLGetRealArray(link, &data,& dims, &heads, &d);
	// 結果を固有値の小さな順にコピー
	// (mathmaticaでは固有値の絶対値の大きな順で出力されるため)
	int i, j, min_j;
	double min_E;
	for(i = 0;i < N;i++){
		min_E = DBL_MAX;
		for(j = 0;j < N;j++){
			if(min_E > ret[j]){
				min_E = ret[j];
				min_j = j;
			};
		};// j
		memcpy(&H[i * N], &data[min_j * N], sizeof(double) * N);
		memcpy(&E[i], &ret[min_j], sizeof(double));
		ret[min_j] = DBL_MAX;
	};
};