求矩阵的秩

       今天编程时需要求矩阵的秩，于是就在网上搜了一个求矩阵的秩的代码(代码源自徐士良《常用算法程序集（C语言描述）》)，用的是全选主元高斯消去法求矩阵的秩。拿了一些实验数据验证了一下，发现有问题。就仔细看了看代码，

 代码如下：

int rank(double* a,int m,int n)  //a为矩阵，m,n分别为矩阵的行与列
{
 int i,j,k,nn,is,js,l,ll,u,v;
    double q,d;
    nn=m;
    if (m>=n) nn=n;
    k=0;
    for (l=0; l<=nn-1; l++)
      { q=0.0;
        for (i=l; i<=m-1; i++)
        for (j=l; j<=n-1; j++)
          { ll=i*n+j; d=fabs(a[ll]);
     if (d>q) { q=d; is=i; js=j;}
          }
        if (q +1.0 == 1.0) return(k);//这里q特别小的时候没有返回，于是我改为if(q< 1e-8)return(k);
        k=k+1;
        if (is!=l)
          { for (j=l; j<=n-1; j++)
              { u=l*n+j; v=is*n+j;
                d=a[u]; a[u]=a[v]; a[v]=d;
              }
          }
        if (js!=l)
          { for (i=l; i<=m-1; i++)
              { u=i*n+js; v=i*n+l;
                d=a[u]; a[u]=a[v]; a[v]=d;
              }
          }
        ll=l*n+l;
        for (i=l+1; i<=n-1; i++)//这里我感觉应该为i<=m-1
          { d=a[i*n+l]/a[ll];
            for (j=l+1; j<=n-1; j++)
              { u=i*n+j;
                a[u]=a[u]-d*a[l*n+j];
              }
          }
      }
    return(k);
  }

不知道改的是否正确，哪位朋友发现改的不对，请指出来，本人不胜感激。