2-2 LU crout 分解

#include <math.h>
#include "..CommonFiles
rutil.h"
#define TINY 1.0e-20;

/* Crout算法
** indx为输出向量，保存部分主元法而改变了行的行排列顺序
** 输出向量d为±1，表示行交换次数为偶数还是为奇数
*/
void ludcmp(float **a, int n, int *indx, float *d) {
    int i,imax,j,k;
    float big,dum,sum,temp;
    float *vv;              /* 用于保存每行的内含比例因子*/

    vv=vector(1,n);

    *d=1.0;
    for (i=1; i<=n; i++) {  /* 按行循环，求内含比例因子(每一行最大元素的倒数)*/
        big=0.0;
        for (j=1; j<=n; j++)
            if ((temp=fabs(a[i][j])) > big) big=temp;
        if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
        vv[i]=1.0/big;
    }

    for (j=1; j<=n; j++) {            /* 外层列循环*/
        for (i=1; i<j; i++) {        /* 上三角部分*/
            sum=a[i][j];
            for (k=1; k<i; k++) 
                sum -= a[i][k]*a[k][j];
            a[i][j]=sum;
        }
        big=0.0;                       /* 初始化最大主元*/
        for (i=j; i<=n; i++) {      /* i >= j 下三角部分*/
            sum=a[i][j];
            for (k=1; k<j; k++)
                sum -= a[i][k]*a[k][j];
            a[i][j]=sum;
            if ( (dum=vv[i]*fabs(sum)) >= big) {
                big=dum;
                imax=i;
            }
        }
        if (j != imax) {            /* 是否进行交换？*/
            for (k=1; k<=n; k++) {
                dum=a[imax][k];
                a[imax][k]=a[j][k];
                a[j][k]=dum;
            }
            *d = -(*d);             /* 改变d的奇偶性*/
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (a[j][j] == 0.0) a[j][j]=TINY;
        if (j != n) {
            dum=1.0/(a[j][j]);
            for (i=j+1; i<=n; i++) a[i][j] *= dum;
        }
    }

    free_vector(vv,1,n);
}
#undef TINY