LSQR - rozwiązanie równań liniowych

Na stronie Stanford University:
http://www.stanford.edu/group/SOL/software/lsqr.html
Znalazłem zaimplementowane rozwiązanie które by mi pasowało gdyby działało.
Nie jestem pewien czy ta funkcja lsqr() nie działa poprawnie, czy to ja nie potrafię jej poprawnie wywołać.
Akurat ten przykład stworzyłem sobie sam, więc wiem jakie jest poprawne rozwiązanie.
Problem polega na tym że funkcja niby znajduje jakieś rozwiązanie zaś z odchyleniem kwadratowym rzędu 1.
A ja chcę znacznie większą dokładność.
W załączniku wszystkie skompresowane źródła.

Podaje swoją próbę odpalenia:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "lsqr.h"

#define Rows 21
#define Cols 9
double A[Rows*Cols]=
  {
   1, 2, 0, 1, 1, 1, 1, 1, 1,
   1, 2, 1, 0, 1, 1, 1, 1, 1,
   1, 2, 1, 1, 1, 1, 0, 1, 1,
   1, 2, 1, 1, 1, 1, 1, 0, 1,
   1, 2, 1, 1, 1, 1, 1, 1, 0,
   1, 1, 2, 0, 1, 1, 1, 1, 1,
   1, 1, 2, 1, 1, 1, 0, 1, 1,
   1, 1, 1, 2, 1, 1, 0, 1, 1,
   1, 0, 1, 1, 2, 1, 1, 1, 1,
   1, 1, 0, 1, 2, 1, 1, 1, 1,
   1, 1, 1, 1, 2, 0, 1, 1, 1,
   1, 1, 1, 1, 2, 1, 0, 1, 1,
   1, 1, 1, 1, 2, 1, 1, 0, 1,
   1, 1, 1, 1, 2, 1, 1, 1, 0,
   1, 0, 1, 1, 1, 2, 1, 1, 1,
   1, 1, 1, 1, 1, 2, 1, 0, 1,
   1, 1, 1, 1, 1, 2, 1, 1, 0,
   1, 1, 1, 1, 1, 1, 1, 2, 0,
   1, 1, 1, 0, 1, 1, 1, 1, 2,
   1, 1, 1, 1, 1, 1, 0, 1, 2,
   1, 1, 1, 0, 2, 1, 1, 1, 1,
  };

double B[Rows]=
  {
    1.6754313690,
    0.5285456570,
    0.6429896060,
    0.5052060890,
    0.4503991730,
   -1.1468855250,
   -1.0324415760,
    0.1144441360,
   -0.4170316480,
    1.2583995340,
    0.0112132130,
    0.2259577710,
    0.0881742540,
    0.0333673380,
   -0.4282446740,
    0.0769612280,
    0.0221543120,
   -0.0548067290,
    0.0781466710,
    0.1925906200,
    0.1115138220,
  };

double GoodX[Cols]=
  {
   -36.31476351,
   5.120326419,
   3.444895237,
   4.591780949,
   4.703294584,
   4.692081558,
   4.477337,
   4.615120517,
   4.669927433
  };

void aprod(int mode,int Row,int Col,double x[],double y[],void *UsrWrk)
  {
   double *A,Sum;
   int c,r,rc,rr;

   A=(double*)UsrWrk;
   if(mode==2)
     {
      for(c=0;c<Col;++c)
        {
         Sum=0;
         for(r=0,rc=c;r<Row;++r,rc+=Col) Sum+=A[rc]*y[r];
         x[c]+=Sum;
        }
     }
   else if(mode==1)
     {
      for(r=0,rr=0;r<Row;++r,rr+=Col)
        {
         Sum=0;
         for(c=0,rc=rr;c<Col;++c,++rc) Sum+=A[rc]*x[c];
         y[r]+=Sum;
        }
     }
  }

int main()
  {
   double *u,*U,*v,*w,*x,*se;
   int istop_out,itn_out,r,c;
   double anorm_out,acond_out,rnorm_out,arnorm_out,xnorm_out,sq,d;

   u=(double*)malloc(Rows*sizeof(double));
   U=(double*)malloc(Rows*sizeof(double));
   v=(double*)malloc(Cols*sizeof(double));
   w=(double*)malloc(Cols*sizeof(double));
   x=(double*)malloc(Cols*sizeof(double));
   se=(double*)malloc(Cols*sizeof(double));
   memcpy(u,B,Rows*sizeof(double));
   memset(U,0,Rows*sizeof(double));
   memset(v,0,Cols*sizeof(double));
   memset(w,0,Cols*sizeof(double));
   memset(x,0,Cols*sizeof(double));
   memset(se,0,Cols*sizeof(double));
   lsqr
     ( 
      Rows,
      Cols,
      aprod,
      0,
      A,
      u,
      v,
      w,
      x,
      se,
      0,
      0,
      0,
      4*Rows*Cols,
      stdout,
      &istop_out,
      &itn_out,
      &anorm_out,
      &acond_out,
      &rnorm_out,
      &arnorm_out,
      &xnorm_out
     );
   printf
    (
     "stop=%d; krokow=%d; acond=%le; anorm=%le; rnorm=%le; xnorm=%le; arnorm=%le;\n\n",
     istop_out,itn_out,acond_out,anorm_out,rnorm_out,xnorm_out,arnorm_out
    );
   printf("%18s %18s\n","ResultRoot","RealRoot");
   for(c=0;c<Cols;++c)
     {
      printf("%18.10f %18.10f\n",x[c],GoodX[c]);
     }
   printf("\n");
   memset(u,0,Cols*sizeof(Rows));
   sq=0;
   aprod(1,Rows,Cols,x,u,A);
   aprod(1,Rows,Cols,GoodX,U,A);
   printf("%18s %18s %18s\n","Result B","Real B","A*RealRoot");
   for(r=0;r<Rows;++r)
     {
      d=u[r]-B[r];
      sq+=d*d;
      printf("%18.10f %18.10f %18.10f\n",u[r],B[r],U[r]);
     }
   printf("Suma kw. %18.10f\n",sq);
  
   getchar();	
   return 0;
  }