c语言编程时出现错误为：1.#QNAN0000000，麻烦大佬们看一下问题出在哪儿了，感谢

Question

数值分析大作业

c语言

这个是数值分析大作业，我是按照书上一步一步编写的，前366行都是没问题的，问题就出在双步位移QR分解那里，但是我真的无能为力了，感谢大佬们相救，感谢。

代码如下：（有些printf是我测试用的）
#include<stdio.h>
#include<math.h>
const int n=10,L=10;
double err=1e-12;
//定义主函数
int main()
{
    int i,j;
    double A[n][n],B[n][n];   //B[][]是为保证矩阵A(n-1)[][]不被破坏的中间矩阵 
    //调用函数声明
    void faketriangle(double a[n][n]);  //构造拟上三角化矩阵函数 
    void QR(double a[n][n]);
    void doublestepQR(double a[n][n]);


    printf("矩阵A为：\n");
    //定义矩阵A

    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {

            if(i==j)
            {
                A[i][j]=1.52*cos(1+i+1.2*(1+j));
            }else
            {
                A[i][j]=sin(0.5*(1+i)+0.2*(1+j));
            }

            printf("%0.12f\t",A[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    printf("矩阵A(n-1)为：\n");

    //调用构造拟上三角化矩阵函数A(n-1)
    faketriangle(A);

    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            B[i][j]=A[i][j];
        }
    }

    printf("\n");
    printf("矩阵A(n-1)分解为QR矩阵\n");
    QR(A);
    printf("\n");
    doublestepQR(B);

 } 



//构造拟上三角化矩阵函数A(n-1) 函数主体
void faketriangle(double a[n][n])
{

    int i,j,r;
    double sum,d,c,h,t;
    double u[n],w[n],q[n],p[n];

    //算法开始迭代 
    for(r=0;r<n-2;r++)
    { 
        sum=0;
        for(i=r+2;i<n;i++)
        {
            sum+=fabs(a[i][r]);
        }

        //判断是否满足a[i][r]=0(i=r+2,..,n)？ 
        if(sum>0)
        {
            sum=0;
            //计算d 
            for(i=r+1;i<n;i++)
            {
                sum+=a[i][r]*a[i][r];
            }
            d=sqrt(sum);

            //计算c 
            if(a[r+1][r]>0)
            {
                c=-d;
            }else
            {
                c=d;
            }
            //计算h
            h=c*c-c*a[r+1][r];


            //向量u[r]的建立 
            for(i=0;i<n;i++)
            {
                if(i<r+1)
                {
                    u[i]=0;
                }else if(i==r+1)
                {
                    u[i]=a[i][r]-c;
                }else
                {
                    u[i]=a[i][r];
                }
            //      printf("%0.12f\n",u[i]);
                }


                //求解向量p 
                for(i=0;i<n;i++)
                {
                    sum=0;
                    for(j=0;j<n;j++)
                    {
                        sum+=a[j][i]*u[j];
                    }
                    p[i]=sum/h;

            //      printf("%0.12f\n",p[i]);
                }
//              
                //求解矩阵q
                for(i=0;i<n;i++)
                {
                    sum=0;
                    for(j=0;j<n;j++)
                    {
                        sum+=a[i][j]*u[j];
                    }
                    q[i]=sum/h;

            //      printf("%0.12f\n",q[i]);
                }

                //求t
                sum=0;
                for(i=0;i<n;i++)
                {
                    sum+=p[i]*u[i];
                }
                t=sum/h;

            //  printf("%0.12f\n",t);

                //求w[]
                for(i=0;i<n;i++)
                {
                    w[i]=q[i]-t*u[i];

            //      printf("%0.12f\n",w[l]);
                }



                //求a(r+1)
                for(i=0;i<n;i++)
                {

                    for(j=0;j<n;j++)
                    {

                        a[i][j]=a[i][j]-(w[i]*u[j]+u[i]*p[j]);

            //          printf("%0.12f\t",a[i][j]);
                    }
            //      printf("\n");
                }
            }
    }

    //输出矩阵A(n-1)
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            printf("%0.12f\t",a[i][j]);
        }
        printf("\n");
    }
} 

//将 A(n-1)矩阵进行QR分解
void QR(double a[n][n])
{
    int i,j,k,r;
    double sum,d,c,h;
    double Q[n][n],u[n],w[n],p[n],R[n][n],b[n][n];

    //定义矩阵Q[][]
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            if(i==j)
            {
                Q[i][j]=1;
            }else
            {
                Q[i][j]=0;
            }
//          printf("%0.12f\t",Q[i][j]);
        }
//      printf("\n");
    }
//  printf("\n");

    //算法迭代
    for(r=0;r<n-1;r++)
    {
        sum=0;
        for(i=r+1;i<n;i++)
        {
            sum+=fabs(a[i][r]);
        }
        if(sum>0)
        {
            //计算d 
            sum=0;
            for(i=r;i<n;i++)
            {
                sum+=a[i][r]*a[i][r];
            }
            d=sqrt(sum);

            //计算c
            if(a[r][r]>0)
            {
                c=-d;
            }else
            {
                c=d;
            }

            //计算h
            h=c*(c-a[r][r]);

            //构造向量u[]
            for(i=0;i<n;i++)
            {
                if(i<r)
                {
                    u[i]=0;
                }else if(i==r)
                {
                    u[i]=a[i][r]-c;
                }else
                {
                    u[i]=a[i][r];
                }
            }

            //计算w[]
            for(i=0;i<n;i++)
            {
                sum=0;
                for(j=0;j<n;j++)
                {
                    sum+=Q[i][j]*u[j];  
                }
                w[i]=sum;
            } 

            //计算Q(r+1)
            for(i=0;i<n;i++)
            {
                for(j=0;j<n;j++)
                {
                    Q[i][j]=Q[i][j]-w[i]*u[j]/h;
                }
            }

            //计算P[]
            for(j=0;j<n;j++)
            {
                sum=0;
                for(i=0;i<n;i++)
                {
                    sum+=a[i][j]*u[i]/h;
                }
                p[j]=sum;
            }

            //计算a[r+1]
            for(i=0;i<n;i++)
            {
                for(j=0;j<n;j++)
                {
                    a[i][j]=a[i][j]-u[i]*p[j];
                }
            }
        }
    }

    printf("Q矩阵为：\n");
    //输出Q[][]、R[][]
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            printf("%0.12f\t",Q[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    printf("R矩阵为：\n");
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            R[i][j]=a[i][j];
            printf("%0.12f\t",R[i][j]);
        }
        printf("\n");
    }

    printf("RQ相乘后的矩阵为：\n");
    //求解R[][]*Q[][]
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            sum=0;
            for(k=0;k<n;k++)
            {
                sum+=R[i][k]*Q[k][j];
            }
            b[i][j]=sum;
            printf("%0.12f\t",b[i][j]);
        }
        printf("\n");
    }

//  printf("QR相乘后的矩阵为：\n");
//  //求解Q[][]*R[][]
//  for(i=0;i<n;i++)
//  {
//      for(j=0;j<n;j++)
//      {
//          sum=0;
//          for(k=0;k<n;k++)
//          {
//              sum+=Q[i][k]*R[k][j];
//          }
//          b[i][j]=sum;
//          printf("%0.12f\t",b[i][j]);
//      }
//      printf("\n");
//  }
}
                                                                       //以上没问题 
//定义复数结构体 
struct complex
{
    double re;
    double im;
};

//对A(n-1)进行双步位移的QR分解
void doublestepQR(double a[n][n])
{
    //M[][]QR分解函数声明 
    void MQR(double a[n][n],double M[n][n],int m);

    int i,j,k,m,r,l;
    double A[n][n],M[n][n],I[n][n],A2[n][n];
    double x,y,z,sum;
    double s,t;

    struct complex lambda[n];
    struct complex s1,s2;

    k=0;m=n-1;

    step3:
    if(fabs(a[m][m-1])<=err)
    {
        lambda[m].re=a[m][m];
        lambda[m].im=0;
        m--;
        goto step4;
    }else
    {
        goto step5;
    }

    step4:
    if(m==0)
    {
        lambda[m].re=a[0][0];
        lambda[m].im=0;
        goto step11;
    }else if(m==-1)
    {
        goto step11;
    }else
    {
        goto step3;
    }

    step5:

    x=a[m-1][m-1]+a[m][m];
    y=a[m-1][m-1]*a[m][m]-a[m-1][m]*a[m][m-1];
    z=x*x-4*y;

    if(z>=0)
    {
        z=sqrt(z); 
        s1.re=(x+z)/2;    
        s1.im=0; 
        s2.re=(x-z)/2;
        s2.im=0;

    }else
    {
        z=sqrt(fabs(z));
        s1.re=(x)/2;
        s1.im=(z)/2;
        s2.re=x/2;
        s2.im=(-z)/2;
    }

    step6:
    if(m==1)
    {
        lambda[m].re=s1.re;
        lambda[m].im=0;
        lambda[m-1].re=s2.re;
        lambda[m-1].im=s2.im;
        goto step11;
    }else
    {
        goto step7;
    }

    step7:
    if(fabs(a[m-1][m-2])<=err)
    {
        if(z>=0)
        {
            lambda[m-1].re=(x+sqrt(z))/2;    //两个特征值
            lambda[m-1].im=0; 
            lambda[m-2].re=(x-sqrt(z))/2;
            lambda[m-2].im=0;

        }else
        {
            lambda[m-1].re=(x)/2;
            lambda[m-1].im=(sqrt(fabs(z)))/2;
            lambda[m-2].re=x/2;
            lambda[m-2].im=(-sqrt(fabs(z)))/2;
        }
        m=m-2;
        goto step4;
    }else
    {
        goto step8;
    }

    step8:
    if(k==L)
    {
        goto step12;

    }else
    {
        goto step9;
    }

    step9:
//  for(i=0;i<m;i++)
//  {
//      for(j=0;j<m;j++)
//      {
//          a[i][j]=a[i][j];
//      }
//  }
    s=a[m-1][m-1]+a[m][m];
    t=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];

    //定义矩阵I[][]
    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            if(i==j)
            {
                I[i][j]=1;
            }else
            {
                I[i][j]=0;
            }
        }
    }
//  printf("%M[][]为：\n"); 
    //计算矩阵M[][]
    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            sum=0;
            for(l=0;l<m;l++)
            {
                sum+=a[i][l]*a[l][j];
            }
            M[i][j]=sum-s*a[i][j]+t*I[i][j];

            printf("%0.12f\t",M[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    //调用M[][]QR分解、计算A(k+1)的函数 
    MQR(a,M,m);

    step10:
        k++;
        goto step3;

    step11:
    printf("特征值已计算完毕。\n");
    for(i=0;i<n;i++)
    {
        printf("%0.12f+%0.12fi\n",lambda[i].re,lambda[i].im);
    }
    for(i=0;i<n;i++)
    {
        if(lambda[i].im==0)
        {
//          vector();
        }
    }

    step12:
    printf("未得到所有特征值。\n"); 

}

//M[][]QR分解、计算A(k+1)的函数 
void MQR(double a[n][n],double M[n][n],int m)
{
    int i,j,r;
    double sum,d,c,h,t;
    double B[n][n],C[n][n],p[n],q[n],v[n],w[n],u[n];

    printf("B:\n");

    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            B[i][j]=M[i][j];
            printf("%f\t",B[i][j]); 
        }
        printf("\n");
    }

    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            C[i][j]=a[i][j];
        }
    }

    //循环计算矩阵A(k+1) 
    for(r=0;r<m;r++)                //r的范围 
    {
        sum=0;
        for(i=r+1;i<m;i++)
        {
            sum+=fabs(B[i][r]);
        }
            printf("sum:\n");
            printf("%0.12f\n",sum);                //sum有问题
            printf("\n\n\n");

        printf("r:%d\n\n\n\n\n",r);

        if(sum>0)
        {
            //计算d
            sum=0;
            for(i=r;i<m+1;i++)
            {
                sum+=B[i][r]*B[i][r];

                printf("B:%0.12f\n",B[i][r]);         //B[][]?
                printf("\n");
            }
            printf("sum:%0.12f\n",sum);                //sum有问题    fiest!!
            printf("\n");

            d=sqrt(sum);
            printf("%0.12f\n",d);                //d有问题
            printf("\n");           
            //计算c
            if(B[r][r]>0)
            {
                c=-d;
            }else
            {
                c=d;
            }
            printf("%0.12f\n",c);              //c有问题 
            printf("\n");           
            //计算h
            h=c*(c-B[r][r]);
            printf("%0.12f\n",h);              //h有问题 
            printf("\n");           
            //构造向量u[]
            for(i=0;i<m;i++)
            {
                if(i<r)
                {
                    u[i]=0;
                }else if(i==r)
                {
                    u[i]=B[i][r]-c;
                }else
                {
                    u[i]=B[i][r];
                }
            }

            //计算v[]
            for(i=0;i<m;i++)
            {
                sum=0;
                for(j=r;j<m;j++)
                {
                    sum+=B[j][i]*u[j];  
                }
                v[i]=sum/h;
            } 

            //计算B(r+1)
            for(i=0;i<m;i++)
            {
                for(j=0;j<m;j++)
                {
                    B[i][j]=B[i][j]-u[i]*v[j];
                }
            }

            //计算p[]
            for(i=0;i<m;i++)
            {
                sum=0;
                for(j=r;j<m;j++)
                {
                    sum+=C[j][i]*u[j];
                }
                p[i]=sum/h;
            }

            //计算q[]
            for(i=0;i<m;i++)
            {
                sum=0;
                for(j=r;j<m;j++)
                {
                    sum+=C[i][j]*u[j];
                }
                q[i]=sum/h;
            }

            //计算t
            sum=0;
            for(i=r;i<m;i++)
            {
                sum+=p[i]*u[i];
            }
            t=sum/h;

            //计算w[]
            for(i=0;i<m;i++)
            {
                w[i]=q[i]-t*u[i];
            } 

            //计算C[r+1]
            for(i=0;i<m;i++)
            {
                for(j=0;j<m;j++)
                {
                    C[i][j]=C[i][j]-w[i]*u[j]-u[i]*p[j];
                }
            }
        }else
        {
            ;
        }
    }

    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            a[i][j]=C[i][j];
            printf("%f\t",a[i][j]);
        }
        printf("\n");
    }
}

Answer 1

for循环出现问题，没有终止条件