对称矩阵的LDLT分解(C语言)

原创 2012年03月30日 18:45:56

/*对称矩阵的LDLT分解
 定理2.2.2:对称矩阵的三角分解:
 设A是n阶对称矩阵,若A的各阶顺序主子式均不等于0,则A可以唯一地分解为
 A=LDL’

-------------A=LDL’的分解算法-------
参考教材:《数值分析》李乃成,梅立泉,科学出版社
    《计算方法教程》第二版 凌永祥,陈明逵

*/

#include<stdio.h>
#include<math.h>

int min(int a,int b);
int main(void)
{
 int size=10;
 int i,j,k,p=0;
 double A[10][10] = {0.0};
 double r[10]={0.0};
 double sum,up=0.0,diag=1.0;
 printf("对称矩阵A[10][10]:\n");
 for(i=0;i<size;i++)
 {
  for(j=0;j<size;j++)
  {
   if(i==j)
    A[i][i]=i+1;
   else
    A[i][j]=min(i,j)-1;
   printf("%5.4f ",A[i][j]);
  }
  printf("\n");
 }
 for(k=0;k<size;k++)
 {
  for(p=0;p<=k-1;p++)
  {
   r[p]=A[p][p]*A[k][p];
  }
  for(p=0;p<=k-1;p++)
  {
   A[k][k]-=A[k][p]*r[p];
  }
  for(i=k+1;i<size;i++)
  {
   sum=0.0;
   for(p=0;p<k;p++)
   {
    sum+=A[i][p]*r[p];
   }
   A[i][k]=(A[i][k]-sum)/A[k][k];
  }
 }

 

 printf("\nLDLT分解后的下三角矩阵L[10][10]:\n");
 for(i=0;i<size;i++)
 {
  for(j=0;j<size;j++)
  {
   if(i<j)
   {
    printf("%5.4f ",up);
   }
   else if(i==j)
   {
    printf("%5.4f ",diag);
   }
   else
   {
    printf("%5.4f ",A[i][j]);
   }
   
  }
  printf("\n");
 }
 printf("\nLDLT分解后的对角阵D[10][10]:\n");
 for(i=0;i<size;i++)
 {
  for(j=0;j<size;j++)
  {
   if(i==j)
   {
    printf("%5.4f ",A[i][i]);
   }
   else
   {
    printf("%5.4f ",up);
   }
   
  }
  printf("\n");
 }

 return 0;
}

int min(int a,int b)
{
 return a<b?a:b;
}

版权声明:本文为博主原创文章,未经博主允许不得转载。

相关文章推荐

LU分解、LDLT分解和Cholesky分解

LU分解 概念:假定我们能把矩阵A写成下列两个矩阵相乘的形式:A=LU,其中L为下三角矩阵,U为上三角矩阵。这样我们可以把线性方程组Ax= b写成 Ax= (LU)x = L(Ux) = b。令Ux ...

QR分解求矩阵特征值、特征向量 C语言

QR分解求矩阵特征值、特征向量            最近在看一个高光谱图像压缩算法,其中涉及到正交变换,计算正交变换时,需要对普通矩阵求其特征向量。想要在网上找一个现成的程序,可能是我百度的能力不强...

[leetcode]Symmetric Tree (对称树 C语言实现)

Symmetric Tree Given a binary tree, check whether it is a mirror of itself (ie, symmetric around i...

c语言的n阶全对称幻方

  • 2011-06-30 20:03
  • 942B
  • 下载

C语言 集合运算 并、交,相对补,对称差,判断两个集合是否相等,求集合幂集

编写程序实现两个集 合的并、交,相对补,对称差,幂集的运算并判断两个集合是否相等
  • Lu_1u
  • Lu_1u
  • 2017-06-27 09:27
  • 666
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:深度学习:神经网络中的前向传播和反向传播算法推导
举报原因:
原因补充:

(最多只允许输入30个字)