-3

I translated some of C code in fortran on my research. While translating I got error in tred2 subroutine. To debug for more easier , I use input for c program while passing in tred2 routine same for fortran. Like (nb =number of band =18 A matrix A(nb,nb) as inout). but I got different output. I use many tred2.f90 subroutine standard and use lapack also for testing. Here is my testing fortran code .. Could you please point out the mistake i made ...

program test_househld
implicit none
INTEGER::nb,i,j
REAL*8:: d1
REAL*8, ALLOCATABLE ::C(:,:),p(:),q(:)

nb =18
ALLOCATE(C(nb,nb),p(nb),q(nb))
do i =1,nb
    p(i) =0.0d0
    q(i) =0.0d0
end do
open(221,file = 'C_matrix.inp',action ='read')
do i =1,nb
        read(221,'(18F24.16)')(C(i,j),j =1,nb)
end do
close(221)

call tred2(C,nb,p,q)

end program test_househld

SUBROUTINE tred2(a,n,d,e)
IMPLICIT NONE
INTEGER :: n
REAL*8 :: a(n,n),d(n),e(n)
INTEGER :: i,j,k,l
REAL*8 ::  f,g,h,hh,scale

DO i=n,2,-1
   l=i-1
   h=0.0D0
   scale=0.0D0
   IF (l > 1) THEN
      scale=SUM(abs(a(i,1:l)))
      IF (scale == 0.0D0) THEN
         e(i)=a(i,l)
      ELSE
         a(i,1:l)=a(i,1:l)/scale
         h=sum(a(i,1:l)**2)
         f=a(i,l)
         g=-sign(sqrt(h),f)
         e(i)=scale*g
         h=h-f*g
         a(i,l)=f-g
         f=0.0D0
         DO j=1,l
            !     Omit following line if finding only eigenvalues
            a(j,i)=a(i,j)/h
            g=0.0D0
            DO k=1,j
               g=g+a(j,k)*a(i,k)
            ENDDO
            DO  k=j+1,l
               g=g+a(k,j)*a(i,k)
            ENDDO
            e(j)=g/h
            f=f+e(j)*a(i,j)
         ENDDO
         hh=f/(h+h)
         DO j=1,l
            f=a(i,j)
            g=e(j)-hh*f
            e(j)=g
            DO k=1,j
               a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
            ENDDO
         ENDDO
      ENDIF
   ELSE
      e(i)=a(i,l)
   ENDIF
   d(i)=h
ENDDO
!     Omit following line if finding only eigenvalues.
d(1)=0.0D0
e(1)=0.0D0
DO i=1,n
   !     Delete lines from here ...
   l=i-1
   IF (d(i) /= 0.0D0) THEN
      DO j=1,l
         g=0.0D0
         DO k=1,l
            g=g+a(i,k)*a(k,j)
         ENDDO
         DO k=1,l
            a(k,j)=a(k,j)-g*a(k,i)
         ENDDO
      ENDDO
   endif
   !     ... to here when finding only eigenvalues.
   d(i)=a(i,i)
   !     Also delete lines from here ...
   a(i,i)=1.
   DO j=1,l
      a(i,j)=0.0D0
      a(j,i)=0.0D0
   ENDDO
   !     ... to here when finding only eigenvalues.
ENDDO

END SUBROUTINE tred2

The used c-code routine is (tred2.c)

#include <math.h>

void tred2(a, n, d, e)
double** a, d[], e[];
int n;
{
    int l, k, j, i;
    double scale, hh, h, g, f;

    for (i = n;i >= 2;i--) {
        l = i - 1;
        h = scale = 0.0;
        if (l > 1) {
            for (k = 1;k <= l;k++)
                scale += fabs(a[i][k]);
            if (scale == 0.0)
                e[i] = a[i][l];
            else {
                for (k = 1;k <= l;k++) {
                    a[i][k] /= scale;
                    h += a[i][k] * a[i][k];
                }
                f = a[i][l];
                g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
                e[i] = scale * g;
                h -= f * g;
                a[i][l] = f - g;
                f = 0.0;
                for (j = 1;j <= l;j++) {
                    a[j][i] = a[i][j] / h;
                    g = 0.0;
                    for (k = 1;k <= j;k++)
                        g += a[j][k] * a[i][k];
                    for (k = j + 1;k <= l;k++)
                        g += a[k][j] * a[i][k];
                    e[j] = g / h;
                    f += e[j] * a[i][j];
                }
                hh = f / (h + h);
                for (j = 1;j <= l;j++) {
                    f = a[i][j];
                    e[j] = g = e[j] - hh * f;
                    for (k = 1;k <= j;k++)
                        a[j][k] -= (f * e[k] + g * a[i][k]);
                }
            }
        }
        else
            e[i] = a[i][l];
        d[i] = h;
    }
    d[1] = 0.0;
    e[1] = 0.0;
    /* Contents of this loop can be omitted if eigenvectors not
                    wanted except for statement d[i]=a[i][i]; */
    for (i = 1;i <= n;i++) {
        l = i - 1;
        if (d[i]) {
            for (j = 1;j <= l;j++) {
                g = 0.0;
                for (k = 1;k <= l;k++)
                    g += a[i][k] * a[k][j];
                for (k = 1;k <= l;k++)
                    a[k][j] -= g * a[k][i];
            }
        }
        d[i] = a[i][i];
        a[i][i] = 1.0;
        for (j = 1;j <= l;j++) a[j][i] = a[i][j] = 0.0;
    }
}

input C_matrix.inp

 -4.9643540965483464      0.0000000002923787     -0.0000000000015816      0.0000000000015866     -0.0000010270946462     -0.0000000000119348     -0.0000014910811312     -0.0000000000083231      0.0000000000051726      0.0000000006747626      0.0000000000042343      0.0000000000044931      0.0000008205316809      0.0000000000013976     -0.0000000000025613      0.0000000000031999      0.0000000000011099     -0.0000002065896463
  0.0000000002923787     -3.0970798965053272      0.0000000000012606     -0.0000000000010828     -0.0000000020105181      0.0000019712291908     -0.0000000019589617     -0.0000000000008715     -0.0000000000035277     -0.0000011866352439     -0.0000000000025535     -0.0000000000030861      0.0000000105614882      0.0000000000007812      0.0000000000003148     -0.0000000000004940     -0.0000000000011315     -0.0000000268806379
 -0.0000000000015816      0.0000000000012606     -3.0726587184630398     -0.0000000000020674     -0.0000000001967980     -0.0000000000031430     -0.0000000053059269      0.0000009717973030     -0.0000000000057488     -0.0000000000018134     -0.0000000000055421     -0.0000000000059327      0.0000000128630816      0.0000000000004427      0.0000000000006127      0.0000005919482648      0.0000000000004205      0.0000000017373891
  0.0000000000015866     -0.0000000000010828     -0.0000000000020674     -3.0726581559742034      0.0000000001970779      0.0000000000026573      0.0000000053540850     -0.0000000000001896      0.0000010540828250      0.0000000000007506      0.0000000000046204      0.0000000000049601     -0.0000000129903302      0.0000000000001229      0.0000000000009025      0.0000000000025886     -0.0000006192561004     -0.0000000019022855
 -0.0000010270946462     -0.0000000020105181     -0.0000000001967980      0.0000000001970779     -1.4799627603758512     -0.0000000000207997      0.0000040022051207      0.0000000000453210      0.0000000000093959     -0.0000000013251917      0.0000000000076454      0.0000000000082252      0.0000012017814372      0.0000000000047051     -0.0000000000118067     -0.0000000000293670     -0.0000000000003220     -0.0000008984357514
 -0.0000000000119348      0.0000019712291908     -0.0000000000031430      0.0000000000026573     -0.0000000000207997     -1.3876784532515820      0.0000000023541494      0.0000000000019625      0.0000000000090133     -0.0000045240623974      0.0000000000124138      0.0000000000078450     -0.0000000018546957      0.0000000000008909     -0.0000000000001635      0.0000000000005688      0.0000000000027158     -0.0000000672281206
 -0.0000014910811312     -0.0000000019589617     -0.0000000053059269      0.0000000053540850      0.0000040022051207      0.0000000023541494     -0.8583134955154931      0.0000000000128576     -0.0000000000192999     -0.0000000014957286     -0.0000000000142530     -0.0000000000152029      0.0000018516481554      0.0000000000039194     -0.0000000000091955     -0.0000000000227099     -0.0000000000003105      0.0000016740054685
 -0.0000000000083231     -0.0000000000008715      0.0000009717973030     -0.0000000000001896      0.0000000000453210      0.0000000000019625      0.0000000000128576     -0.6260557686606371      0.0000000000034069     -0.0000000000008690      0.0000000000023244      0.0000000000025902     -0.0000000015017165     -0.0000000000038208      0.0000000000039194      0.0000046551165327      0.0000000000227413     -0.0000000010756788
  0.0000000000051726     -0.0000000000035277     -0.0000000000057488      0.0000010540828250      0.0000000000093959      0.0000000000090133     -0.0000000000192999      0.0000000000034069     -0.6260536476274911      0.0000000000115542      0.0000000000137463      0.0000000000147947     -0.0000000024681981      0.0000000000006085     -0.0000000000004013      0.0000000000248150     -0.0000047190361606     -0.0000000002257375
  0.0000000006747626     -0.0000011866352439     -0.0000000000018134      0.0000000000007506     -0.0000000013251917     -0.0000045240623974     -0.0000000014957286     -0.0000000000008690      0.0000000000115542     -0.5803541891855615      0.0000000000049024      0.0000000000096723     -0.0000000023349175      0.0000000000037884      0.0000000000025921     -0.0000000000039961      0.0000000000017223      0.0000004168417969
  0.0000000000042343     -0.0000000000025535     -0.0000000000055421      0.0000000000046204      0.0000000000076454      0.0000000000124138     -0.0000000000142530      0.0000000000023244      0.0000000000137463      0.0000000000049024     -0.5574486988431754      0.0000000000119368     -0.0000000012816454      0.0000063881906259      0.0000000008948734      0.0000000000000565      0.0000000000038825     -0.0000000002143914
  0.0000000000044931     -0.0000000000030861     -0.0000000000059327      0.0000000000049601      0.0000000000082252      0.0000000000078450     -0.0000000000152029      0.0000000000025902      0.0000000000147947      0.0000000000096723      0.0000000000119368     -0.5574485377065610     -0.0000000013589895     -0.0000000008438915      0.0000063940897453      0.0000000000005915      0.0000000000039688     -0.0000000002034879
  0.0000008205316809      0.0000000105614882      0.0000000128630816     -0.0000000129903302      0.0000012017814372     -0.0000000018546957      0.0000018516481554     -0.0000000015017165     -0.0000000024681981     -0.0000000023349175     -0.0000000012816454     -0.0000000013589895     -0.4561210997669952      0.0000000000031373     -0.0000000000086804     -0.0000000000561001     -0.0000000000078430     -0.0000006303634354
  0.0000000000013976      0.0000000000007812      0.0000000000004427      0.0000000000001229      0.0000000000047051      0.0000000000008909      0.0000000000039194     -0.0000000000038208      0.0000000000006085      0.0000000000037884      0.0000063881906259     -0.0000000008438915      0.0000000000031373     -0.2436725624222180     -0.0000000000051047     -0.0000000000123756      0.0000000000002341     -0.0000000011331551
 -0.0000000000025613      0.0000000000003148      0.0000000000006127      0.0000000000009025     -0.0000000000118067     -0.0000000000001635     -0.0000000000091955      0.0000000000039194     -0.0000000000004013      0.0000000000025921      0.0000000008948734      0.0000063940897453     -0.0000000000086804     -0.0000000000051047     -0.2436723950073047      0.0000000000031071      0.0000000000049463      0.0000000018149829
  0.0000000000031999     -0.0000000000004940      0.0000005919482648      0.0000000000025886     -0.0000000000293670      0.0000000000005688     -0.0000000000227099      0.0000046551165327      0.0000000000248150     -0.0000000000039961      0.0000000000000565      0.0000000000005915     -0.0000000000561001     -0.0000000000123756      0.0000000000031071     -0.2191906816051025     -0.0000000000013001     -0.0000000003356968
  0.0000000000011099     -0.0000000000011315      0.0000000000004205     -0.0000006192561004     -0.0000000000003220      0.0000000000027158     -0.0000000000003105      0.0000000000227413     -0.0000047190361606      0.0000000000017223      0.0000000000038825      0.0000000000039688     -0.0000000000078430      0.0000000000002341      0.0000000000049463     -0.0000000000013001     -0.2191900460492452      0.0000000002003037
 -0.0000002065896463     -0.0000000268806379      0.0000000017373891     -0.0000000019022855     -0.0000008984357514     -0.0000000672281206      0.0000016740054685     -0.0000000010756788     -0.0000000002257375      0.0000004168417969     -0.0000000002143914     -0.0000000002034879     -0.0000006303634354     -0.0000000011331551      0.0000000018149829     -0.0000000003356968      0.0000000002003037     -0.0602903370568763
Raj Kumar Paudel
  • 121
  • 3
  • 2
    Sorry, it looks like you're trying to get us to debug your program. That's not what SO is for. What's more, you just claim that you get 'different output', but don't tell us what it should be, what it is, and how it is different. And you're also not showing us the c program. All of these things make it impossible for us to help you even if we wanted to. – chw21 Dec 01 '20 at 04:14
  • Hi @chw21, yeah i added the c -code here. Actually I supplied same input for both c and fortran program but got different result. I try to understand numerical recipe book tred2.f, they used physical dimension Np, but in my case the matrix is 18*18 band so i put N and Np equal to nb. I am looking for help to solve this issue .. – Raj Kumar Paudel Dec 01 '20 at 04:31
  • Edit the question to provide a [mre] including full output from both the FORTRAN and C versions with enough precision to uniquely distinguish all numbers (at least 17 significant digits). – Eric Postpischil Dec 01 '20 at 09:04
  • Isn't it not so strange that the result is affected rather strongly by round-off errors (if the matrix elements differ so much)? I've tried comparing the above code (from numerical recipes?), tred2.f [here](http://www.netlib.no/netlib/eispack/tred2.f), and some julia's routine, and the results seem to agree up to ~8 significant digits... – roygvib Dec 01 '20 at 12:49

1 Answers1

0

There is no reason to expect the results to be the same. The error is in your expectation. As just one reason the output could be different (and there are many), compare:

        for (k = 1;k <= j;k++)
            a[j][k] -= (f * e[k] + g * a[i][k]);

And

        DO k=1,j
           a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
        ENDDO

In floating point math, a-(b+c) is not the same as (a-b)-c. This is particularly the case if b+c or a-b is different in magnitude from a or c. Where precision is lost affects the results.

Here's some code to demonstrate one way it can fail:

#include <stdio.h>

void compare(double a, double b, double c)
{
   double d = a-(b+c);
   double e = (a-b)-c;
   printf("d=%lf e=%lf d-e=%lf d-e=0->%s\n",
       d, e, d-e, ((d-e) == 0) ? "yes" : "no");
}

int main()
{
    compare (1.0, 2.0, 3.0);                    // happens to be the same
    compare (1e18L, 1e18L-0.0001L, 0.00001);    // happens to not be the same
}

See the output here.

David Schwartz
  • 179,497
  • 17
  • 214
  • 278
  • ,I change the expectation bracket accordingly but the result not different. After 8 some place after 10 decimal point other the result differs i tried to debug not able to solve this issue... one more i only read upto 15 digit and then i use QP to read 16 digit but the output has after 8 or 10 decimal point differnt... – Raj Kumar Paudel Dec 01 '20 at 05:22
  • @RajKumarPaudel Did you look at the example? One result was zero and the other was -0.00001 -- that's a difference in five decimal places. Unless you can show that the result is *wrong*, you don't have anything unexpected. Different in expected. – David Schwartz Dec 01 '20 at 05:27
  • I checked it on the other hand my fortran code and C code are both form numerical recipe on C and fortran 77... I guess the problem is dimension. In numerical recipe fortran ..SUBROUTINE tred2(a,n,np,d,e) (the input n and Np and tell Np >=n) in my case i use same n=Np =nb. I do not know where is the problem.. whie reading input 15 digit read correctly but I use QP to get 16 digit...but not able to solve this issue – Raj Kumar Paudel Dec 01 '20 at 05:41