[Mplapack-devel] Bug in Rgerfs

[Fletcher, J. P <j.p...@as...>]

// solve_dd_refine.cpp from  solve_double_refine.cpp
// solve_dd.cpp from inv_mat_dd.cpp
//
// Get inverse of Matrix A via Rgetri, using double
// This file is freely usable.
// written by Nakata Maho, 2009/9/23.

#include <mblas_dd.h>
#include <mlapack_dd.h>
#include <stdio.h>

//Matlab/Octave format
void printmat(int N, int M, dd_real *A, int LDA)
{
    dd_real mtmp;

    printf("[ ");
    for (int i = 0; i < N; i++) {
	printf("[ ");
	for (int j = 0; j < M; j++) {
	    mtmp = A[i + j * LDA];
	    std::cout << mtmp;
	    //	    printf("%5.2e", mtmp);
	    if (j < M - 1)
		printf(", ");
	}
	if (i < N - 1)
	    printf("]; ");
	else
	    printf("] ");
    }
    printf("]");
}

int main()
{
    mpackint n = 3;
    mpackint lwork, info;

    dd_real *A = new dd_real[n * n];
    dd_real *AF = new dd_real[n * n];
    mpackint *ipiv = new mpackint[n];
    mpackint *ipiv1 = new mpackint[n+1];
    mpackint *iwork = new mpackint[n+1];

    mpackint nrhs = 1;
 
    dd_real *B = new dd_real[n*nrhs];
    dd_real *X = new dd_real[n*nrhs];
    dd_real *B1 = new dd_real[n*nrhs+1];
    dd_real *X1 = new dd_real[n*nrhs+1];
    dd_real *work = new dd_real[3 * n + 1];
    dd_real *berr = new dd_real[n];
    dd_real *ferr = new dd_real[n];
//setting A matrix
    A[0 + 0 * n] = 1;    A[0 + 1 * n] = 4;    A[0 + 2 * n] = 6;
    A[1 + 0 * n] = 2;    A[1 + 1 * n] = 9;    A[1 + 2 * n] = 14;
    A[2 + 0 * n] = 3;    A[2 + 1 * n] = 14;   A[2 + 2 * n] = 23;
 
    B[0] = 1; B[1] = 2; B[2] = 3;

    printf("A = ");
    printmat(n, n, A, n);
    printf("\n");
    printf("B = ");
    printmat(n, nrhs, B, n);
    printf("\n");

    for (int i = 0; i < n*n; i++) AF[i] = A[i];
    for (int i = 0; i < n*nrhs; i++) X[i] = B[i];

    Rgetrf(n, n, AF, n, ipiv, &info);
    Rgetrs("n", n, nrhs, AF, n, ipiv, X, n, &info);

    printf("X = ");
    printmat(n, nrhs, X, n);
    printf("\n");
 
    for (int i = 0; i < n; i++) ipiv1[i+1] = ipiv[i];
    for (int i = 0; i < n; i++) B1[i+1] = B[i];
    for (int i = 0; i < n; i++) X1[i+1] = X[i];

    Rgerfs("n", n, nrhs, A, n, AF, n, ipiv1, B1, n, X1, n,
           ferr, berr, work, iwork, &info);
    printf("info = %d\n",info);
    printf("revised X = ");
    printmat(n, nrhs, X, n);
    printf("\n");
    printmat(nrhs, 1, ferr, n);
    printmat(nrhs, 1, berr, n);
    printf("\n");

    delete[]work;
    delete[]iwork;
    delete[]ipiv;
    delete[]ipiv1;
    delete[]A;
    delete[]AF;
    delete[]B;
    delete[]X;
    delete[]berr;
    delete[]ferr;


}