Cholesky分解法可以将矩阵分解为,其中L为_半正定矩阵cholesky分解

Cholesky分解法可以将矩阵分解为,其中L为_半正定矩阵cholesky分解头文件:/**Copyright(c)2008-2011ZhangMing(M.Zhang),zmjerry@163.com**Thisprogramisfreesoftware;youcanredistributeitand/ormodifyit*underthetermsoftheGNUGeneralPublicLicenseasp…

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全家桶1年46,售后保障稳定

头文件:

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com

*

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

*

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

*

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

*

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

*/

/*****************************************************************************

* cholesky.h

*

* Class template of Cholesky decomposition.

*

* For a symmetric, positive definite matrix A, this function computes the

* Cholesky factorization, i.e. it computes a lower triangular matrix L such

* that A = L*L’. If the matrix is not symmetric or positive definite, the

* function computes only a partial decomposition. This can be tested with

* the isSpd() flag.

*

* This class also supports factorization of complex matrix by specializing

* some member functions.

*

* Adapted form Template Numerical Toolkit.

*

* Zhang Ming, 2010-01 (revised 2010-12), Xi’an Jiaotong University.

*****************************************************************************/

#ifndef CHOLESKY_H

#define CHOLESKY_H

#include

namespace splab

{

template

class Cholesky

{

public:

Cholesky();

~Cholesky();

bool isSpd() const;

void dec( const Matrix &A );

Matrix getL() const;

Vector solve( const Vector &b );

Matrix solve( const Matrix &B );

private:

bool spd;

Matrix L;

};

//class Cholesky

#include

}

// namespace splab

#endif

// CHOLESKY_H

实现文件:

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com

*

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

*

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

*

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

*

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

*/

/*****************************************************************************

* cholesky-impl.h

*

* Implementation for Cholesky class.

*

* Zhang Ming, 2010-01 (revised 2010-12), Xi’an Jiaotong University.

*****************************************************************************/

/**

* constructor and destructor

*/

template

Cholesky::Cholesky() : spd(true)

{

}

template

Cholesky::~Cholesky()

{

}

/**

* return true, if original matrix is symmetric positive-definite.

*/

template

inline bool Cholesky::isSpd() const

{

return spd;

}

/**

* Constructs a lower triangular matrix L, such that L*L’= A.

* If A is not symmetric positive-definite (SPD), only a partial

* factorization is performed. If isspd() evalutate true then

* the factorizaiton was successful.

*/

template

void Cholesky::dec( const Matrix &A )

{

int m = A.rows();

int n = A.cols();

spd = (m == n);

if( !spd )

return;

L = Matrix(n,n);

// main loop.

for( int j=0; j

{

Type d = 0;

for( int k=0; k

{

Type s = 0;

for( int i=0; i

s += L[k][i]*L[j][i];

L[j][k] = s = (A[j][k]-s) / L[k][k];

d = d + s*s;

spd = spd && (A[k][j] == A[j][k]);

}

d = A[j][j] – d;

spd = spd && ( d > 0 );

L[j][j] = sqrt( d > 0 ? d : 0 );

for( int k=j+1; k

L[j][k] = 0;

}

}

/**

* return the lower triangular factor, L, such that L*L’=A.

*/

template

inline Matrix Cholesky::getL() const

{

return L;

}

/**

* Solve a linear system A*x = b, using the previously computed

* cholesky factorization of A: L*L’.

*/

template

Vector Cholesky::solve( const Vector &b )

{

int n = L.rows();

if( b.dim() != n )

return Vector();

Vector x = b;

// solve L*y = b

for( int k=0; k

{

for( int i=0; i

x[k] -= x[i]*L[k][i];

x[k] /= L[k][k];

}

// solve L’*x = y

for( int k=n-1; k>=0; –k )

{

for( int i=k+1; i

x[k] -= x[i]*L[i][k];

x[k] /= L[k][k];

}

return x;

}

/**

* Solve a linear system A*X = B, using the previously computed

* cholesky factorization of A: L*L’.

*/

template

Matrix Cholesky::solve( const Matrix &B )

{

int n = L.rows();

if( B.rows() != n )

return Matrix();

Matrix X = B;

int nx = B.cols();

// solve L*Y = B

for( int j=0; j

for( int k=0; k

{

for( int i=0; i

X[k][j] -= X[i][j]*L[k][i];

X[k][j] /= L[k][k];

}

// solve L’*X = Y

for( int j=0; j

for( int k=n-1; k>=0; –k )

{

for( int i=k+1; i

X[k][j] -= X[i][j]*L[i][k];

X[k][j] /= L[k][k];

}

return X;

}

/**

* Main loop of specialized member function. This macro definition is

* aimed at avoiding code duplication.

*/

#define MAINLOOP \

{ \

int m = A.rows(); \

int n = A.cols(); \

\

spd = (m == n); \

if( !spd ) \

return; \

\

for( int j=0; j

{ \

spd = spd && (imag(A[j][j]) == 0); \

d = 0; \

\

for( int k=0; k

{ \

s = 0; \

for( int i=0; i

s += L[k][i] * conj(L[j][i]); \

\

L[j][k] = s = (A[j][k]-s) / L[k][k]; \

d = d + norm(s); \

spd = spd && (A[k][j] == conj(A[j][k])); \

} \

\

d = real(A[j][j]) – d; \

spd = spd && ( d > 0 ); \

\

L[j][j] = sqrt( d > 0 ? d : 0 ); \

for( int k=j+1; k

L[j][k] = 0; \

} \

}

/**

* Solving process of specialized member function. This macro definition is

* aimed at avoiding code duplication.

*/

#define SOLVE1 \

{ \

for( int k=0; k

{ \

for( int i=0; i

x[k] -= x[i]*L[k][i]; \

\

x[k] /= L[k][k]; \

} \

\

for( int k=n-1; k>=0; –k ) \

{ \

for( int i=k+1; i

x[k] -= x[i]*conj(L[i][k]); \

\

x[k] /= L[k][k]; \

} \

\

return x; \

}

/**

* Solving process of specialized member function. This macro definition is

* aimed at avoiding code duplication.

*/

#define SOLVE2 \

{ \

int nx = B.cols(); \

for( int j=0; j

for( int k=0; k

{ \

for( int i=0; i

X[k][j] -= X[i][j]*L[k][i]; \

\

X[k][j] /= L[k][k]; \

} \

\

for( int j=0; j

for( int k=n-1; k>=0; –k ) \

{ \

for( int i=k+1; i

X[k][j] -= X[i][j]*conj(L[i][k]); \

\

X[k][j] /= L[k][k]; \

} \

\

return X; \

}

/**

* Specializing for “dec” member function.

*/

template <>

void Cholesky >::dec( const Matrix > &A )

{

float d;

complex s;

L = Matrix >(A.cols(),A.cols());

MAINLOOP;

}

template <>

void Cholesky >::dec( const Matrix > &A )

{

double d;

complex s;

L = Matrix >(A.cols(),A.cols());

MAINLOOP;

}

template <>

void Cholesky >::dec( const Matrix > &A )

{

long double d;

complex s;

L = Matrix >(A.cols(),A.cols());

MAINLOOP;

}

/**

* Specializing for “solve” member function.

*/

template <>

Vector > Cholesky >::solve( const Vector > &b )

{

int n = L.rows();

if( b.dim() != n )

return Vector >();

Vector > x = b;

SOLVE1

}

template <>

Vector > Cholesky >::solve( const Vector > &b )

{

int n = L.rows();

if( b.dim() != n )

return Vector >();

Vector > x = b;

SOLVE1

}

template <>

Vector > Cholesky >::solve( const Vector > &b )

{

int n = L.rows();

if( b.dim() != n )

return Vector >();

Vector > x = b;

SOLVE1

}

/**

* Specializing for “solve” member function.

*/

template <>

Matrix > Cholesky >::solve( const Matrix > &B )

{

int n = L.rows();

if( B.rows() != n )

return Matrix >();

Matrix > X = B;

SOLVE2

}

template <>

Matrix > Cholesky >::solve( const Matrix > &B )

{

int n = L.rows();

if( B.rows() != n )

return Matrix >();

Matrix > X = B;

SOLVE2

}

template <>

Matrix > Cholesky >::solve( const Matrix > &B )

{

int n = L.rows();

if( B.rows() != n )

return Matrix >();

Matrix > X = B;

SOLVE2

}

测试代码:

/*****************************************************************************

* cholesky_test.cpp

*

* Cholesky class testing.

*

* Zhang Ming, 2010-01 (revised 2010-12), Xi’an Jiaotong University.

*****************************************************************************/

#define BOUNDS_CHECK

#include

#include

#include

using namespace std;

using namespace splab;

typedef double Type;

const int N = 5;

int main()

{

Matrix A(N,N), L(N,N);

Vector b(N);

for( int i=1; i

{

for( int j=1; j

if( i == j )

A(i,i) = i;

else

if( i < j )

A(i,j) = i;

else

A(i,j) = j;

b(i) = i*(i+1)/2.0 + i*(N-i);

}

cout << setiosflags(ios::fixed) << setprecision(3);

cout << “The original matrix A : ” << A << endl;

Cholesky cho;

cho.dec(A);

if( !cho.isSpd() )

cout << “Factorization was not complete.” << endl;

else

{

L = cho.getL();

cout << “The lower triangular matrix L is : ” << L << endl;

cout << “A – L*L^T is : ” << A – L*trT(L) << endl;

Vector x = cho.solve(b);

cout << “The constant vector b : ” << b << endl;

cout << “The solution of Ax = b : ” << x << endl;

cout << “The Ax – b : ” << A*x-b << endl;

Matrix IA = cho.solve(eye(N,Type(1)));

cout << “The inverse matrix of A : ” << IA << endl;

cout << “The product of A*inv(A) : ” << A*IA << endl;

}

return 0;

}

运行结果:

The original matrix A : size: 5 by 5

1.000 1.000 1.000 1.000 1.000

1.000 2.000 2.000 2.000 2.000

1.000 2.000 3.000 3.000 3.000

1.000 2.000 3.000 4.000 4.000

1.000 2.000 3.000 4.000 5.000

The lower triangular matrix L is : size: 5 by 5

1.000 0.000 0.000 0.000 0.000

1.000 1.000 0.000 0.000 0.000

1.000 1.000 1.000 0.000 0.000

1.000 1.000 1.000 1.000 0.000

1.000 1.000 1.000 1.000 1.000

A – L*L^T is : size: 5 by 5

0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000

0.000 0.000 0.000 0.000 0.000

The constant vector b : size: 5 by 1

5.000

9.000

12.000

14.000

15.000

The solution of Ax = b : size: 5 by 1

1.000

1.000

1.000

1.000

1.000

The Ax – b : size: 5 by 1

0.000

0.000

0.000

0.000

0.000

The inverse matrix of A : size: 5 by 5

2.000 -1.000 0.000 0.000 0.000

-1.000 2.000 -1.000 0.000 0.000

0.000 -1.000 2.000 -1.000 0.000

0.000 0.000 -1.000 2.000 -1.000

0.000 0.000 0.000 -1.000 1.000

The product of A*inv(A) : size: 5 by 5

1.000 0.000 0.000 0.000 0.000

0.000 1.000 0.000 0.000 0.000

0.000 0.000 1.000 0.000 0.000

0.000 0.000 0.000 1.000 0.000

0.000 0.000 0.000 0.000 1.000

Process returned 0 (0x0) execution time : 0.109 s

Press any key to continue.

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/219123.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

  • java中hashcode的用法_javahashcode作用

    java中hashcode的用法_javahashcode作用hashcode()是干什么用的?首先hashcode是哈希算法的一中简单实现,他是一个对象的哈希吗值。一般和equals一起使用。 hashcode也是用来查找的,如果你学过数据结构就应该知道,在查找和排序这一章有 例如内存中有这样的位置 01234567 而我有个类,这个类有个字段叫ID,我要把这个类存放在以上8个位置之一,如果不用hashcode而任意存放,

  • Vue之Axios跨域问题解决方案

    Vue之Axios跨域问题解决方案背景:因为axios中只能使用get和post方法来进行请求数据,没有提供jsonp等方法进行跨域访问数据axios中文网址:https://www.kancloud.cn/yunye/axios/234845//axios中的GET请求axios.get(‘/user’,{params:{ID:‘001’}})…

  • Ti杯电子竞赛前期准备工作

    Ti杯电子竞赛前期准备工作标题竞赛时间和竞赛周期:1997年开始每二年举办-届,竞赛时间定于竞赛举办年度的9月份,赛期四天三夜(具体日期届时通知)。在双数的非竞赛年份,组织开展全国的专题性竞赛。竞赛方式:全国统一-命题、分赛区组织的方式。采用“半封闭、相对集中”的组织方式进行。学生可查阅纸介或网络技术资料,队内学生集体商讨设计,分工负责、团结协作,以队为基本单位独立完成竞赛任务;竞赛期间不允许任何教师或其他…

  • shuffle洗牌算法java_洗牌算法shuffle[通俗易懂]

    洗牌算法1.背景阿里的面试的时候做的一道笔试题:题目:写一个方法,入参为自然数n(n>0),返回一个自然数数组,数组长度为n,元素为[1,n]之间,且每个元素不重复,数组中各元素顺序要求随机;实例1:输入:N=3输出:132实例2:输入:N=5输出:32514当时我的解法(写了两种方法):写的好烂,面完和面试官交流的时候面试官让我看下Collect…

  • php fastcgi,配置apache以fastcgi运行php[通俗易懂]

    php fastcgi,配置apache以fastcgi运行php[通俗易懂]apache默认是用自带的mod_php模块运行php,现在我们介绍使用fastcgi来执行php脚本。先说下fastcgi的优点:Fastcgi的优点:从稳定性上看,fastcgi是以独立的进程池运行来cgi,单独一个进程死掉,系统可以很轻易的丢弃,然后重新分配新的进程来运行逻辑.·从安全性上看,Fastcgi支持分布式运算.fastcgi和宿主的server完全独立,fastc…

  • Python学习系列之lambda表达式

    Python学习系列之lambda表达式一、lambda定义与用法lambda表达式是一行的函数。它们在其他语言中也被称为匿名函数。即,函数没有具体的名称,而用def创建的方法是有名称的。如果你不想在程序中对一个函数使用两次,你也许会想用lambda表达式,它们和普通的函数完全一样。而且当使用函数作为参数的时候,lambda表达式非常有用,可以让代码简单,简洁。lambda表达式返回的是function类型,说明是一个函数类型。…

    2022年10月18日

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号