//奇异值分解法求广义逆
//本函数返回值小于0表示在奇异值分解过程,
//中迭代值超过了60次还未满足精度要求.
//返回值大于0表示正常返回。
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//aa-长度为n*m的数组，返回式存放A的广义逆
//eps-精度要求
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav()
int dginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka)
{ int i,j,k,l,t,p,q,f;
    i=dluav(a,m,n,u,v,eps,ka);
    if (i<0)
{
  return(-1);
}
    j=n;
    if (m<n)
{
  j=m;
}
    j=j-1;
    k=0;
    while ((k<=j)&&(a[k*n+k]!=0.0)) 
{
  k=k+1;
}
    k=k-1;
    for (i=0; i<=n-1; i++)
{
  for (j=0; j<=m-1; j++)
  { 
   t=i*m+j; 
   aa[t]=0.0;
   for (l=0; l<=k; l++)
   { 
    f=l*n+i; 
    p=j*m+l; 
    q=l*n+l;
    aa[t]=aa[t]+v[f]*u[p]/a[q];
   }
  }
}
    return(1);
} //实数矩阵的奇异值分解