改写UMFPACK算例中的压缩方式（二）

两次遍历，第一次遍历可以确定非零元素的个数；第二次遍历确定位置和存储非零元

//Data:2013-2-24
//修改了Ai Ax的类型 利用最大维数n*n来保存，可以调用正确结果 不过不知系统随机分配的值 函数没有用
//Data:2013-2-26
//调用UMFPACK包来实现求解方程组
//UMFPACK采用CSC（列压缩存储） matlab中的接口为A/b
#include <stdio.h>
#include <math.h>
#include "umfpack.h"
//传递的四个参数A b x n -----Data:2013-02-27
//意思为Ax=b n为矩阵维数
int umf(double **A,double *b,double *x,int n)
{
    //printf("--");
    //printf("%d \n",A[1][1]);
    // int n=5;
    double *null =(double *)NULL ;
    void *Symbolic, *Numeric ;
    int i,j;

/*    
    //定义矩阵A
    int **A=new int *[n];
    for(j=0;j<n;j++)
        A[j]= new int[n];

    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            scanf("%d",&A[i][j]);
*/

    int *Ap=new int [n+1];
    Ap[0]=0;

   /* int *Ai=new int[n*n];
    for(i=0;i<n*n;i++)
        Ai[i]=1;
    double *Ax=new double[n*n];
    for(i=0;i<n*n;i++)
        Ax[i]=1;*/
    double epsilon=0.00001;
    int NZnum=0;//矩阵非零元的个数

    int k=0,l=0,m=0;
    for(j=0;j<n;j++)
    {
        for(i=0;i<n;i++)
        {
            if(abs(A[i][j])>epsilon)
            {
             /*   Ai[l++]=i;
                Ax[m++]=A[i][j];*/


                NZnum++;
            }
        }
        Ap[k+1]=NZnum;
        k++;
    }
    Ap[n]=NZnum;

     int *Ai=new int[NZnum];
    double *Ax=new double[NZnum];

    // int k=0,l=0,m=0;
    for(j=0;j<n;j++)
    {
        for(i=0;i<n;i++)
        {
            if(abs(A[i][j])>epsilon)
            {
                Ai[l++]=i;
                Ax[m++]=A[i][j];
            }
        }
    }

    //printf("--Ai---\n");
    //for (int i=0; i<n*n; i++)
    //{
    //    printf("%d ", Ai[i]);
    //}
    //printf("\n");
    
    umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null);
    umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null) ;
    umfpack_di_free_symbolic (&Symbolic);
    umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null);
    umfpack_di_free_numeric (&Numeric) ;

    //for(i=0;i<n;i++) 
    //    printf("x[%d]=%g\n", i, x[i]) ;
    return 0;
}