文档章节

线性代数方程组的求解

RapidBird
 RapidBird
发布于 2010/03/30 13:02
字数 4047
阅读 440
收藏 1

#include "stdafx.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "LinearEquation.h"
#include "MatrixAlgo.h"

//求解三对角线方程组的追赶法
int atrde(double b[],int n, int m, double d[])
{
 int k,j;
 double s;
 if(m != (3*n-2))
 {
  printf("err\n");
  return(-2);
 }
 for(k = 0;k <= n-2;k++)
 {
  j = 3*k;
  s = b[j];
  if(fabs(s)+1.0 == 1.0)
  {
   printf("fail\n");
   return(0);
  }
  b[j+1] = b[j+1]/s;
  d[k] = d[k]/s;
  b[j+3] = b[j+3]-b[j+2]*b[j+1];
  d[k+1] = d[k+1]-b[j+2]*d[k];
 }
 s = b[3*n-3];
 if(fabs(s)+1.0 == 1.0)
 {
  printf("fail\n");
  return(0);
 }
 d[n-1] = d[n-1]/s;
 for(k = n-2;k>=0;k--)
  d[k] = d[k]-b[3*k+1]*d[k+1];
 return(2);
}

//求解托伯利兹方程组的列文逊方法
int atlvs(double t[], int n, double b[], double x[])
{
 int i,j,k;
 double a,beta,q,c,h,*y,*s;
 s = (double *)malloc(n*sizeof(double));
 y = (double *)malloc(n*sizeof(double));
 a = t[0];
 if(fabs(a)+1.0 == 1.0)
 {
  free(s);
  free(y);
  printf("fail\n");
  return(-1);
 }
 y[0] = 1.0;
 x[0] = b[0]/a;
 for(k = 1; k <= n-1; k++)
 {
  beta = 0.0;
  q = 0.0;
  for(j = 0; j <= k-1; j++)
  {
   beta = beta+y[j]*t[j+1];
   q = q+x[j]*t[k-j];
  }
  if(fabs(a)+1.0 == 1.0)
  {
   free(s);
   free(y);
   printf("fail\n");
   return(-1);
  }
  c = -beta/a;
  s[0] = c*y[k-1];
  y[k] = y[k-1];
  if(k != 1)
   for(i = 1; i <= k-1; i++)
    s[i] = y[i-1]+c*y[k-i-1];
  a = a+c*beta;
  if(fabs(a)+1.0 == 1.0)
  {
   free(s);
   free(y);
   printf("fail\n");
   return(-1);
  }
  h = (b[k]-q)/a;
  for(i = 0; i <= k-1; i++)
  {
   x[i] = x[i]+h*s[i];
   y[i] = s[i];
  }
  x[k] = h*y[k];
 }
 free(s);
 free(y);
 return(1);
}

//求解对称方程组的分解法
int aldle(double a[], int n, int m, double c[])
{
 int i,j,l,k,u,v,w,k1,k2,k3;
 double p;
 if(fabs(a[0])+1.0 == 1.0)
 {
  printf("fail\n");
  return(-2);
 }
 for(i = 1; i <= n-1; i++)
 {
  u = i*n;
  a[u] = a[u]/a[0];
 }
 for(i = 1; i <= n-2; i++)
 {
  u = i*n+i;
  for(j = 1; j <= i; j++)
  {
   v = i*n+j-1;
   l = (j-1)*n+j-1;
   a[u] = a[u]-a[v]*a[v]*a[l];
  }
  p = a[u];
  if(fabs(p)+1.0 == 1.0)
  {
   printf("fail\n");
   return(-2);
  }
  for(k = i+1; k <= n-1; k++)
  {
   u = k*n+i;
   for(j = 1; j <= i; j++)
   {
    v = k*n+j-1;
    l = i*n+j-1;
    w = (j-1)*n+j-1;
    a[u] = a[u]-a[v]*a[l]*a[w];
   }
   a[u] = a[u]/p;
  }
 }
 u = n*n-1;
 for(j = 1; j <= n-1; j++)
 {
  v = (n-1)*n+j-1;
  w = (j-1)*n+j-1;
  a[u] = a[u]-a[v]*a[v]*a[w];
 }
 p = a[u];
 if(fabs(p)+1.0 == 1.0)
 {
  printf("fail\n");
  return(-2);
 }
 for(j = 0; j <= m-1; j++)
  for(i = 1; i <= n-1; i++)
  {
   u = i*m+j;
   for(k = 1; k <= i; k++)
   {
    v = i*n+k-1;
    w = (k-1)*m+j;
    c[u] = c[u]-a[v]*c[w];
   }
  }
  for(i = 1; i <= n-1; i++)
  {
   u = (i-1)*n+i-1;
   for(j = i; j <= n-1; j++)
   {
    v = (i-1)*n+j;
    w = j*n+i-1;
    a[v] = a[u]*a[w];
   }
  }
  for(j = 0; j <= m-1; j++)
  {
   u = (n-1)*m+j;
   c[u] = c[u]/p;
   for(k = 1; k <= n-1; k++)
   {
    k1 = n-k;
    k3 = k1-1;
    u = k3*m+j;
    for(k2 = k1; k2 <= n-1; k2++)
    {
     v = k3*n+k2;
     w = k2*m+j;
     c[u] = c[u]-a[v]*c[w];
    }
    c[u] = c[u]/a[k3*n+k3];
   }
  }
 return(2);
}

//高斯-赛德尔迭代法
int agsdl(double a[], double b[], int n, double x[], double eps)
{
 int i,j,u,v;
 double p,t,s,q;
 for(i = 0; i <= n-1; i++)
 {
  u = i*n+i;
  p = 0.0;
  x[i] = 0.0;
  for(j = 0; j <= n-1; j++)
   if(i != j)
   {
    v = i*n+j;
    p = p+fabs(a[v]);
   }
  if(p>=fabs(a[u]))
  {
   printf("fail\n");
   return(-1);
  }
 }
 p = eps+1.0;
 while(p>=eps)
 {
  p = 0.0;
  for(i = 0; i <= n-1; i++)
  {
   t = x[i];
   s = 0.0;
   for(j = 0; j <= n-1; j++)
    if(j != i)
     s = s+a[i*n+j]*x[j];
   x[i] = (b[i]-s)/a[i*n+i];
   q = fabs(x[i]-t)/(1.0+fabs(x[i]));
   if(q>p)
    p = q;
  }
 }
 return(1);
}

//求解对称正定方程组的共轭梯度法
void agrad(double a[], int n, double b[], double eps, double x[])
{
 int i,k;
 double *p,*r,*s,*q,alpha,beta,d,e;
 p = (double *)malloc(n*sizeof(double));
 r = (double *)malloc(n*sizeof(double));
 s = (double *)malloc(n*sizeof(double));
 q = (double *)malloc(n*sizeof(double));
 for(i = 0; i <= n-1; i++)
 {
  x[i] = 0.0;
  p[i] = b[i];
  r[i] = b[i];
 }
 i = 0;
 while(i <= n-1)
 {
  brmul(a,p,n,n,1,s);
  d = 0.0;
  e = 0.0;
  for(k = 0; k <= n-1; k++)
  {
   d = d+p[k]*b[k];
   e = e+p[k]*s[k];
  }
  alpha = d/e;
  for(k = 0; k <= n-1; k++)
   x[k] = x[k]+alpha*p[k];
  brmul(a,x,n,n,1,q);
  d = 0.0;
  for(k = 0; k <= n-1; k++)
  {
   r[k] = b[k]-q[k];
   d = d+r[k]*s[k];
  }
  beta = d/e;
  d = 0.0;
  for(k = 0; k <= n-1; k++)
   d = d+r[k]*r[k];
  d = sqrt(d);
  if(d<eps)
  {
   free(p);
   free(r);
   free(s);
   free(q);
   return;
  }
  for(k = 0; k <= n-1; k++)
   p[k] = r[k]-beta*p[k];
  i = i+1;
 }
 free(p);
 free(r);
 free(s);
 free(q);
 return;
}

//求解线性最小二乘问题的豪斯荷尔德变换法
int agmqr(double a[], int m, int n, double b[], double q[])
{
 int i,j;
 double d,*c;
 c = (double *)malloc(n*sizeof(double));
 i = bmaqr(a,m,n,q);
 if(i == 0)
 {
  free(c);
  return(0);
 }
 for(i = 0; i <= n-1; i++)
 {
  d = 0.0;
  for(j = 0; j <= m-1; j++)
   d = d+q[j*m+i]*b[j];
  c[i] = d;
 }
 b[n-1] = c[n-1]/a[n*n-1];
 for(i = n-2; i>=0; i--)
 {
  d = 0.0;
  for(j = i+1; j <= n-1; j++)
   d = d+a[i*n+j]*b[j];
  b[i] = (c[i]-d)/a[i*n+i];
 }
 free(c);
 return(1);
}

//求解线性最小二乘问题的广义逆法
int agmiv(double a[], int m, int n, double b[], double x[], double aa[], double eps, double u[], double v[], int ka)
{
 int i,j;
 i = bginv(a,m,n,aa,eps,u,v,ka);
 if(i<0)
  return(-1);
 for(i = 0; i <= n-1; i++)
 {
  x[i] = 0.0;
  for(j = 0; j <= m-1; j++)
   x[i] = x[i]+aa[i*m+j]*b[j];
 }
 return(1);
}

//全选主元高斯-约当消去法
int agjdn(double a[], double b[], int n, int m)
{
 int *js,l,k,i,j,is,p,q;
 double d,t;
 js = (int *)malloc(n*sizeof(int));
 l = 1;
 for(k = 0;k <= n-1;k++)
 {
  d = 0.0;
  for(i = k;i <= n-1;i++)
   for(j = k;j <= n-1;j++)
   {
    t = fabs(a[i*n+j]);
    if(t>d)
    {
     d = t;
     js[k] = j;
     is = i;
    }
   }
  if(d+1.0 == 1.0)
   l = 0;
  else
  {
   if(js[k] != k)
    for(i = 0;i <= n-1;i++)
    {
     p = i*n+k;
     q = i*n+js[k];
     t = a[p];
     a[p] = a[q];
     a[q] = t;
    }
   if(is != k)
   {
    for(j = k;j <= n-1;j++)
    {
     p = k*n+j;
     q = is*n+j;
     t = a[p];
     a[p] = a[q];
     a[q] = t;
    }
    for(j = 0;j <= m-1;j++)
    {
     p = k*m+j;
     q = is*m+j;
     t = b[p];
     b[p] = b[q];
     b[q] = t;
    }
   }
  }
  if(l == 0)
  {
   free(js);
   printf("fail\n");
   return(0);
  }
  d = a[k*n+k];
  for(j = k+1;j <= n-1;j++)
  {
   p = k*n+j;
   a[p] = a[p]/d;
  }
  for(j = 0;j <= m-1;j++)
  {
   p = k*m+j;
   b[p] = b[p]/d;
  }
  for(j = k+1;j <= n-1;j++)
   for(i = 0;i <= n-1;i++)
   {
    p = i*n+j;
    if(i != k)
     a[p] = a[p]-a[i*n+k]*a[k*n+j];
   }
  for(j = 0;j <= m-1;j++)
   for(i = 0;i <= n-1;i++)
   {
    p = i*m+j;
    if(i != k)
     b[p] = b[p]-a[i*n+k]*b[k*m+j];
   }
 }
 for(k = n-1;k>=0;k--)
  if(js[k] != k)
   for(j = 0;j <= m-1;j++)
   {
    p = k*m+j;
    q = js[k]*m+j;
    t = b[p];
    b[p] = b[q];
    b[q] = t;
   }
 free(js);
 return(1);
}

//求解大型稀疏方程组的全选主元高斯-约当消去法
int aggje(double a[], int n, double b[])
{
 int *js,i,j,k,is,u,v;
 double d,t;
 js = (int *)malloc(n*sizeof(int));
 for(k = 0; k <= n-1; k++)
 {
  d = 0.0;
  for(i = k; i <= n-1; i++)
   for(j = k; j <= n-1; j++)
   {
    t = fabs(a[i*n+j]);
    if(t>d)
    {
     d = t;
     js[k] = j;
     is = i;
    }
   }
  if(d+1.0 == 1.0)
  {
   free(js);
   printf("fail\n");
   return(0);
  }
  if(is != k)
  {
   for(j = k; j <= n-1; j++)
   {
    u = k*n+j;
    v = is*n+j;
    t = a[u];
    a[u] = a[v];
    a[v] = t;
   }
   t = b[k];
   b[k] = b[is];
   b[is] = t;
  }
  if(js[k] != k)
   for(i=0; i <= n-1; i++)
   {
    u = i*n+k;
    v = i*n+js[k];
    t = a[u];
    a[u] = a[v];
    a[v] = t;
   }
  t = a[k*n+k];
  for(j = k+1; j <= n-1; j++)
  {
   u = k*n+j;
   if(a[u] != 0.0)
    a[u] = a[u]/t;
  }
  b[k] = b[k]/t;
  for(j = k+1; j <= n-1; j++)
  {
   u = k*n+j;
   if(a[u] != 0.0)
   {
    for(i = 0; i <= n-1; i++)
    {
     v = i*n+k;
     if((i != k)&&(a[v] != 0.0))
     {
      is = i*n+j;
      a[is] = a[is]-a[v]*a[u];
     }
    }
   }
  }
  for(i = 0; i <= n-1; i++)
  {
   u = i*n+k;
   if((i != k)&&(a[u] != 0.0))
    b[i] = b[i]-a[u]*b[k];
  }
 }
 for(k = n-1; k>=0; k--)
  if(k != js[k])
   {
    t = b[k];
    b[k] = b[js[k]];
    b[js[k]] = t;
   }
 free(js);
 return(1);
}

//全选主元高斯消去法
int agaus(double a[], double b[], int n)
{
 int *js,l,k,i,j,is,p,q;
 double d,t;
 js = (int *)malloc(n*sizeof(int));
 l = 1;
 for(k = 0;k <= n-2;k++)
 {
  d = 0.0;
  for(i = k;i <= n-1;i++)
   for(j = k;j <= n-1;j++)
   {
    t = fabs(a[i*n+j]);
    if(t>d)
    {
     d = t;
     js[k] = j;
     is = i;
    }
   }
  if(d+1.0 == 1.0)
   l = 0;
  else
  {
   if(js[k] != k)
   for(i = 0;i <= n-1;i++)
   {
    p = i*n+k;
    q = i*n+js[k];
    t = a[p];
    a[p] = a[q];
    a[q] = t;
   }
   if(is != k)
   {
    for(j = k;j <= n-1;j++)
    {
     p = k*n+j;
     q = is*n+j;
     t = a[p];
     a[p] = a[q];
     a[q] = t;
    }
    t = b[k];
    b[k] = b[is];
    b[is] = t;
   }
  }
  if(l == 0)
  {
   free(js);
   printf("fail\n");
   return(0);
  }
  d = a[k*n+k];
  for(j = k+1;j <= n-1;j++)
  {
   p = k*n+j;
   a[p] = a[p]/d;
  }
  b[k] = b[k]/d;
  for(i = k+1;i <= n-1;i++)
  {
   for(j = k+1;j <= n-1;j++)
   {
    p = i*n+j;
    a[p] = a[p]-a[i*n+k]*a[k*n+j];
   }
   b[i] = b[i]-a[i*n+k]*b[k];
  }
 }
 d = a[(n-1)*n+n-1];
 if(fabs(d)+1.0 == 1.0)
 {
  free(js);
  printf("fail\n");
  return(0);
 }
 b[n-1] = b[n-1]/d;
 for(i = n-2;i>=0;i--)
 {
  t = 0.0;
  for(j = i+1;j <= n-1;j++)
   t = t+a[i*n+j]*b[j];
  b[i] = b[i]-t;
 }
 js[n-1] = n-1;
 for(k = n-1;k>=0;k--)
  if(js[k] != k)
  {
   t = b[k];
   b[k] = b[js[k]];
   b[js[k]] = t;
  }
 free(js);
 return(1);
}

//复系数方程组的全选主元高斯-约当消去法
int acjdn(double ar[], double ai[], double br[], double bi[], int n, int m)
{
 int *js,l,k,i,j,is,u,v;
 double p,q,s,d;
 js = (int *)malloc(n*sizeof(int));
 for(k = 0;k <= n-1;k++)
 {
  d = 0.0;
  for(i = k;i <= n-1;i++)
   for(j = k;j <= n-1;j++)
   {
    u = i*n+j;
    p = ar[u]*ar[u]+ai[u]*ai[u];
    if(p>d)
    {
     d = p;
     js[k] = j;
     is = i;
    }
   }
  if(d+1.0 == 1.0)
  {
   free(js);
   printf("err**fail\n");
   return(0);
  }
  if(is != k)
  {
   for(j = k;j <= n-1;j++)
   {
    u = k*n+j;
    v = is*n+j;
    p = ar[u];
    ar[u] = ar[v];
    ar[v] = p;
    p = ai[u];
    ai[u] = ai[v];
    ai[v] = p;
   }
   for(j = 0;j <= m-1;j++)
   {
    u = k*m+j;
    v = is*m+j;
    p = br[u];
    br[u] = br[v];
    br[v] = p;
    p = bi[u];
    bi[u] = bi[v];
    bi[v] = p;
   }
  }
  if(js[k] != k)
   for(i = 0;i <= n-1;i++)
   {
    u = i*n+k;
    v = i*n+js[k];
    p = ar[u];
    ar[u] = ar[v];
    ar[v] = p;
    p = ai[u];
    ai[u] = ai[v];
    ai[v] = p;
   }
  v = k*n+k;
  for(j = k+1;j <= n-1;j++)
  {
   u = k*n+j;
   p = ar[u]*ar[v];
   q = -ai[u]*ai[v];
   s = (ar[v]-ai[v])*(ar[u]+ai[u]);
   ar[u] = (p-q)/d;
   ai[u] = (s-p-q)/d;
  }
  for(j = 0;j <= m-1;j++)
  {
   u = k*m+j;
   p = br[u]*ar[v];
   q = -bi[u]*ai[v];
   s = (ar[v]-ai[v])*(br[u]+bi[u]);
   br[u] = (p-q)/d;
   bi[u] = (s-p-q)/d;
  }
  for(i = 0;i <= n-1;i++)
   if(i != k)
   {
    u = i*n+k;
    for(j = k+1;j <= n-1;j++)
    {
     v = k*n+j;
     l = i*n+j;
     p = ar[u]*ar[v];
     q = ai[u]*ai[v];
     s = (ar[u]+ai[u])*(ar[v]+ai[v]);
     ar[l] = ar[l]-p+q;
     ai[l] = ai[l]-s+p+q;
    }
    for(j = 0;j <= m-1;j++)
    {
     l = i*m+j;
     v = k*m+j;
     p = ar[u]*br[v];
     q = ai[u]*bi[v];
     s = (ar[u]+ai[u])*(br[v]+bi[v]);
     br[l] = br[l]-p+q;
     bi[l] = bi[l]-s+p+q;
    }
   }
 }
 for(k = n-1;k>=0;k--)
  if(js[k] != k)
   for(j = 0;j <= m-1;j++)
   {
    u = k*m+j;
    v = js[k]*m+j;
    p = br[u];
    br[u] = br[v];
    br[v] = p;
    p = bi[u];
    bi[u] = bi[v];
    bi[v] = p;
   }
 free(js);
 return(1);
}

//求解对称正定方程组的平方根法
int achol(double a[], int n, int m, double d[])
{
 int i,j,k,u,v;
 if((a[0]+1.0 == 1.0)||(a[0]<0.0))
 {
  printf("fail\n");
  return(-2);
 }
 a[0] = sqrt(a[0]);
 for(j = 1; j <= n-1; j++)
  a[j] = a[j]/a[0];
 for(i = 1; i <= n-1; i++)
 {
  u = i*n+i;
  for(j = 1; j <= i; j++)
  {
   v = (j-1)*n+i;
   a[u] = a[u]-a[v]*a[v];
  }
  if((a[u]+1.0 == 1.0)||(a[u]<0.0))
  {
   printf("fail\n");
   return(-2);
  }
  a[u] = sqrt(a[u]);
  if(i != (n-1))
  {
   for(j = i+1; j <= n-1; j++)
   {
    v = i*n+j;
    for(k = 1; k <= i; k++)
     a[v] = a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
    a[v] = a[v]/a[u];
   }
  }
 }
 for(j = 0; j <= m-1; j++)
 {
  d[j] = d[j]/a[0];
  for(i = 1; i <= n-1; i++)
  {
   u = i*n+i;
   v = i*m+j;
   for(k = 1; k <= i; k++)
    d[v] = d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
   d[v] = d[v]/a[u];
  }
 }
 for(j = 0; j <= m-1; j++)
 {
  u = (n-1)*m+j;
  d[u] = d[u]/a[n*n-1];
  for(k = n-1; k>=1; k--)
  {
   u = (k-1)*m+j;
   for(i = k; i <= n-1; i++)
   {
    v = (k-1)*n+i;
    d[u] = d[u]-a[v]*d[i*m+j];
   }
   v = (k-1)*n+k-1;
   d[u] = d[u]/a[v];
  }
 }
 return(2);
}

//复系数方程组的全选主元高斯消去法
int acgas(double ar[], double ai[], int n, double br[], double bi[])
{
 int *js,l,k,i,j,is,u,v;
 double p,q,s,d;
 js = (int *)malloc(n*sizeof(int));
 for(k = 0;k <= n-2;k++)
 {
  d = 0.0;
  for(i = k;i <= n-1;i++)
   for(j = k;j <= n-1;j++)
   {
    u = i*n+j;
    p = ar[u]*ar[u]+ai[u]*ai[u];
    if(p>d)
    {
     d = p;
     js[k] = j;
     is = i;
    }
   }
  if(d+1.0 == 1.0)
  {
   free(js);
   printf("err**fail\n");
   return(0);
  }
  if(is != k)
  {
   for(j = k;j <= n-1;j++)
   {
    u = k*n+j;
    v = is*n+j;
    p = ar[u];
    ar[u] = ar[v];
    ar[v] = p;
    p = ai[u];
    ai[u] = ai[v];
    ai[v] = p;
   }
   p = br[k];
   br[k] = br[is];
   br[is] = p;
   p = bi[k];
   bi[k] = bi[is];
   bi[is] = p;
  }
  if(js[k] != k)
   for(i = 0;i <= n-1;i++)
   {
    u = i*n+k;
    v = i*n+js[k];
    p = ar[u];
    ar[u] = ar[v];
    ar[v] = p;
    p = ai[u];
    ai[u] = ai[v];
    ai[v] = p;
   }
  v = k*n+k;
  for(j = k+1;j <= n-1;j++)
  {
   u = k*n+j;
   p = ar[u]*ar[v];
   q = -ai[u]*ai[v];
   s = (ar[v]-ai[v])*(ar[u]+ai[u]);
   ar[u] = (p-q)/d;
   ai[u] = (s-p-q)/d;
  }
  p = br[k]*ar[v];
  q = -bi[k]*ai[v];
  s = (ar[v]-ai[v])*(br[k]+bi[k]);
  br[k] = (p-q)/d;
  bi[k] = (s-p-q)/d;
  for(i = k+1;i <= n-1;i++)
  {
   u = i*n+k;
   for(j = k+1;j <= n-1;j++)
   {
    v = k*n+j;
    l = i*n+j;
    p = ar[u]*ar[v];
    q = ai[u]*ai[v];
    s = (ar[u]+ai[u])*(ar[v]+ai[v]);
    ar[l] = ar[l]-p+q;
    ai[l] = ai[l]-s+p+q;
   }
   p = ar[u]*br[k];
   q = ai[u]*bi[k];
   s = (ar[u]+ai[u])*(br[k]+bi[k]);
   br[i] = br[i]-p+q;
   bi[i] = bi[i]-s+p+q;
  }
 }
 u = (n-1)*n+n-1;
 d = ar[u]*ar[u]+ai[u]*ai[u];
 if(d+1.0 == 1.0)
 {
  free(js);
  printf("err**fail\n");
  return(0);
 }
 p = ar[u]*br[n-1];
 q = -ai[u]*bi[n-1];
 s = (ar[u]-ai[u])*(br[n-1]+bi[n-1]);
 br[n-1] = (p-q)/d;
 bi[n-1] = (s-p-q)/d;
 for(i = n-2;i>=0;i--)
  for(j = i+1;j <= n-1;j++)
  {
   u = i*n+j;
   p = ar[u]*br[j];
   q = ai[u]*bi[j];
   s = (ar[u]+ai[u])*(br[j]+bi[j]);
   br[i] = br[i]-p+q;
   bi[i] = bi[i]-s+p+q;
  }
 js[n-1] = n-1;
 for(k = n-1;k>=0;k--)
  if(js[k] != k)
  {
   p = br[k];
   br[k] = br[js[k]];
   br[js[k]] = p;
   p = bi[k];
   bi[k] = bi[js[k]];
   bi[js[k]] = p;
  }
 free(js);
 return(1);
}

//病态方程组的求解
int abint(double a[], int n, double b[], double eps, double x[])
{
 int i,j,k,kk;
 double *p,*r,*e,q,qq;
 p = (double *)malloc(n*n*sizeof(double));
 r = (double *)malloc(n*sizeof(double));
 e = (double *)malloc(n*sizeof(double));
 i = 60;
 for(k = 0; k <= n-1; k++)
  for(j = 0; j <= n-1; j++)
   p[k*n+j] = a[k*n+j];
 for(k = 0; k <= n-1; k++)
  x[k] = b[k];
 kk = agaus(p,x,n);
 if(kk == 0)
 {
  free(p);
  free(r);
  free(e);
  return(kk);
 }
 q = 1.0+eps;
 while(q>=eps)
 {
  if(i == 0)
  {
   free(p);
   free(r);
   free(e);
   return(i);
  }
  i = i-1;
  brmul(a,x,n,n,1,e);
  for( k = 0; k <= n-1; k++)
   r[k] = b[k]-e[k];
  for( k = 0; k <= n-1; k++)
   for( j = 0; j <= n-1; j++)
    p[k*n+j] = a[k*n+j];
  kk = agaus(p,r,n);
  if(kk == 0)
  {
   free(p);
   free(r);
   free(e);
   return(kk);
  }
  q = 0.0;
  for( k = 0; k <= n-1; k++)
  {
   qq = fabs(r[k])/(1.0+fabs(x[k]+r[k]));
   if(qq>q)
    q = qq;
  }
  for( k = 0; k <= n-1; k++)
   x[k] = x[k]+r[k];
 }
 free(p);
 free(r);
 free(e);
 return(1);
}

//一般带型方程组的求解
int aband(double b[], double d[], int n, int l, int il, int m)
{
 int ls,k,i,j,is,u,v;
 double p,t;
 if(il != (2*l+1))
 {
  printf("fail\n");
  return(-2);
 }
 ls = l;
 for(k = 0;k <= n-2;k++)
 {
  p = 0.0;
  for(i = k;i <= ls;i++)
  {
   t = fabs(b[i*il]);
   if(t>p)
   {
    p = t;
    is = i;
   }
  }
  if(p+1.0 == 1.0)
  {
   printf("fail\n");
   return(0);
  }
  for(j = 0;j <= m-1;j++)
  {
   u = k*m+j;
   v = is*m+j;
   t = d[u];
   d[u] = d[v];
   d[v] = t;
  }
  for(j = 0;j <= il-1;j++)
  {
   u = k*il+j;
   v = is*il+j;
   t = b[u];
   b[u] = b[v];
   b[v] = t;
  }
  for(j = 0;j <= m-1;j++)
  {
   u = k*m+j;
   d[u] = d[u]/b[k*il];
  }
  for(j = 1;j <= il-1;j++)
  {
   u = k*il+j;
   b[u] = b[u]/b[k*il];
  }
  for(i = k+1;i <= ls;i++)
  {
   t = b[i*il];
   for(j = 0;j <= m-1;j++)
   {
    u = i*m+j;
    v = k*m+j;
    d[u] = d[u]-t*d[v];
   }
   for(j = 1;j <= il-1;j++)
   {
    u = i*il+j;
    v = k*il+j;
    b[u-1] = b[u]-t*b[v];
   }
   u = i*il+il-1;
   b[u] = 0.0;
  }
  if(ls != (n-1))
   ls = ls+1;
 }
 p = b[(n-1)*il];
 if(fabs(p)+1.0 == 1.0)
 {
  printf("fail\n");
  return(0);
 }
 for(j = 0;j <= m-1;j++)
 {
  u = (n-1)*m+j;
  d[u] = d[u]/p;
 }
 ls = 1;
 for(i = n-2;i>=0;i--)
 {
  for(k = 0;k <= m-1;k++)
  {
   u = i*m+k;
   for(j = 1;j <= ls;j++)
   {
    v = i*il+j;
    is = (i+j)*m+k;
    d[u] = d[u]-b[v]*d[is];
   }
  }
  if(ls != (il-1))
   ls = ls+1;
 }
 return(2);
}

                                     ----根据《C语言常用算法程序集》整理

<投票>

 

© 著作权归作者所有

上一篇: 计算机程序设计
下一篇: 特殊函数
RapidBird
粉丝 11
博文 24
码字总数 64196
作品 0
海淀
私信 提问
有限元分析软件--OpenFEM

有限元分析,即使用有限元方法来分析静态或动态的物体或系统。在这种方法中一个物体或系统被分解为由多个相互联结的、简单、独立的点组成的几何模型。在这 种方法中这些独立的点的数量是有限...

匿名
2009/10/27
11.6K
0
[DEEP LEARNING An MIT Press book in preparation]Linear algebra

线性代数是数学的一个重要分支,经常被应用到工程问题中,要理解深度学习以及操作深度学习,那么对于线性代数深刻的理解是非常重要的,以下摘要是我从DL book的第二章线性代数中抽取出来的比...

sunbaigui
2014/10/10
0
0
使用javascript实现矩阵LU分解

版权声明:本文为博主原创文章,欢迎读者转载学习,转载请注明文章出处。 https://blog.csdn.net/qq_37338983/article/details/79918254 在线性代数中,LU分解是将一个矩阵分解为 L(单位下三...

Bruce_wjh
2018/04/12
0
0
Function Set in OPEN CASCADE

Function Set in OPEN CASCADE eryar@163.com Abstract. The common math algorithms library provides a C++ implementation of the most frequently used mathematical algorithms. These ......

eryar
2016/01/16
84
0
如何使用python进行常规方程求解(Sympy or Scipy)(非线性方程组,一元二次方程,多元一次方程,因式分解等)

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/meiqi0538/article/details/82990432 前言 在科学计算中,我们经常会遇到数值计算,可能遇到高数,线性代数等...

皮乾东
2018/10/09
0
0

没有更多内容

加载失败,请刷新页面

加载更多

Giraph源码分析(八)—— 统计每个SuperStep中参与计算的顶点数目

作者|白松 目的:科研中,需要分析在每次迭代过程中参与计算的顶点数目,来进一步优化系统。比如,在SSSP的compute()方法最后一行,都会把当前顶点voteToHalt,即变为InActive状态。所以每次...

数澜科技
50分钟前
1
0
Xss过滤器(Java)

问题 最近旧的系统,遇到Xss安全问题。这个系统采用用的是spring mvc的maven工程。 解决 maven依赖配置 <properties><easapi.version>2.2.0.0</easapi.version></properties><dependenci......

亚林瓜子
今天
7
0
Navicat 快捷键

操作 结果 ctrl+q 打开查询窗口 ctrl+/ 注释sql语句 ctrl+shift +/ 解除注释 ctrl+r 运行查询窗口的sql语句 ctrl+shift+r 只运行选中的sql语句 F6 打开一个mysql命令行窗口 ctrl+l 删除一行 ...

低至一折起
今天
7
0
Set 和 Map

Set 1:基本概念 类数组对象, 内部元素唯一 let set = new Set([1, 2, 3, 2, 1]); console.log(set); // Set(3){ 1, 2, 3 } [...set]; // [1, 2, 3] 接收数组或迭代器对象 ...

凌兮洛
今天
1
0
PyTorch入门笔记一

张量 引入pytorch,生成一个随机的5x3张量 >>> from __future__ import print_function>>> import torch>>> x = torch.rand(5, 3)>>> print(x)tensor([[0.5555, 0.7301, 0.5655],......

仪山湖
今天
5
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部