Thread: [Mplapack-devel] Bug in Rgerfs
Status: Pre-Alpha
Brought to you by:
nakatamaho
From: Fletcher, J. P <j.p...@as...> - 2011-04-06 08:18:57
|
// 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; } |
From: Maho N. <ch...@ma...> - 2011-04-07 00:35:04
|
Hi Fletcher, committed, thanks! Best regards, Nakata Maho From: "Fletcher, John P" <j.p...@as...> Subject: Bug in Rgerfs Date: Wed, 06 Apr 2011 09:18:41 +0100 > Nakata Maho > > I have several things to send you and will start with this one. > > I have been testing Rgerfs for double and dd and did not get the results I get with LAPACK. > > I think the reason is that some of the accesses in Rgerfs to ipiv, B and X start at [1] and not [0], and also access to work is wrongly indexed so that one more variable is accessed at the end. > > I attach two test routines. solve_double_refine.cpp shows it not working and solve_dd_refine.cpp shows reindexing the ipiv, B and X arrays - I copied them to new arrays - and extra work. This works except that the error results are wrong - file attached. > > The other file is the correct results for the same problem using the routine from LAPACK. > > I hope this helps. > > John Fletcher > > P.S. I will send other things later. > > Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678 > Chemical Engineering and Applied Chemistry (CEAC), > Associate Dean - External Relations, > School of Engineering and Applied Science (EAS), > Aston University, Aston Triangle, BIRMINGHAM B4 7ET U.K. > > |