开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 数学计算 > 查看源码
Matrix.cpp
资源名称：科学与工程数值计算算法配套源码vc++.zip　[点击查看]
上传用户：weigute
上传日期：2007-03-02
资源大小：1287k
文件大小：71k
源码类别：

数学计算

开发平台：

Visual C++

Matrix.cpp：源码内容
  1. //////////////////////////////////////////////////////////////////////
  2. // Matrix.cpp
  3. //
  4. // 操作矩阵的类 CMatrix 的实现文件
  5. //
  6. // 周长发编制, 2002/8
  7. //////////////////////////////////////////////////////////////////////
  9. #include "stdafx.h"
  10. #include "Matrix.h"
  12. #ifdef _DEBUG
  13. #undef THIS_FILE
  14. static char THIS_FILE[]=__FILE__;
  15. #define new DEBUG_NEW
  16. #endif
  18. //////////////////////////////////////////////////////////////////////
  19. // Construction/Destruction
  20. //////////////////////////////////////////////////////////////////////
  22. //////////////////////////////////////////////////////////////////////
  23. // 基本构造函数
  24. //////////////////////////////////////////////////////////////////////
  25. CMatrix::CMatrix()
  26. {
  27. m_nNumColumns = 1;
  28. m_nNumRows = 1;
  29. m_pData = NULL;
  30. BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
  31. ASSERT(bSuccess);
  32. }
  34. //////////////////////////////////////////////////////////////////////
  35. // 指定行列构造函数
  36. //
  37. // 参数：
  38. // 1. int nRows - 指定的矩阵行数
  39. // 2. int nCols - 指定的矩阵列数
  40. //////////////////////////////////////////////////////////////////////
  41. CMatrix::CMatrix(int nRows, int nCols)
  42. {
  43. m_nNumRows = nRows;
  44. m_nNumColumns = nCols;
  45. m_pData = NULL;
  46. BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
  47. ASSERT(bSuccess);
  48. }
  50. //////////////////////////////////////////////////////////////////////
  51. // 指定值构造函数
  52. //
  53. // 参数：
  54. // 1. int nRows - 指定的矩阵行数
  55. // 2. int nCols - 指定的矩阵列数
  56. // 3. double value[] - 一维数组，长度为nRows*nCols，存储矩阵各元素的值
  57. //////////////////////////////////////////////////////////////////////
  58. CMatrix::CMatrix(int nRows, int nCols, double value[])
  59. {
  60. m_nNumRows = nRows;
  61. m_nNumColumns = nCols;
  62. m_pData = NULL;
  63. BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
  64. ASSERT(bSuccess);
  66. SetData(value);
  67. }
  69. //////////////////////////////////////////////////////////////////////
  70. // 方阵构造函数
  71. //
  72. // 参数：
  73. // 1. int nSize - 方阵行列数
  74. //////////////////////////////////////////////////////////////////////
  75. CMatrix::CMatrix(int nSize)
  76. {
  77. m_nNumRows = nSize;
  78. m_nNumColumns = nSize;
  79. m_pData = NULL;
  80. BOOL bSuccess = Init(nSize, nSize);
  81. ASSERT (bSuccess);
  82. }
  84. //////////////////////////////////////////////////////////////////////
  85. // 方阵构造函数
  86. //
  87. // 参数：
  88. // 1. int nSize - 方阵行列数
  89. // 2. double value[] - 一维数组，长度为nRows*nRows，存储方阵各元素的值
  90. //////////////////////////////////////////////////////////////////////
  91. CMatrix::CMatrix(int nSize, double value[])
  92. {
  93. m_nNumRows = nSize;
  94. m_nNumColumns = nSize;
  95. m_pData = NULL;
  96. BOOL bSuccess = Init(nSize, nSize);
  97. ASSERT (bSuccess);
  99. SetData(value);
  100. }
  102. //////////////////////////////////////////////////////////////////////
  103. // 拷贝构造函数
  104. //
  105. // 参数：
  106. // 1. const CMatrix& other - 源矩阵
  107. //////////////////////////////////////////////////////////////////////
  108. CMatrix::CMatrix(const CMatrix& other)
  109. {
  110. m_nNumColumns = other.GetNumColumns();
  111. m_nNumRows = other.GetNumRows();
  112. m_pData = NULL;
  113. BOOL bSuccess = Init(m_nNumRows, m_nNumColumns);
  114. ASSERT(bSuccess);
  116. // copy the pointer
  117. memcpy(m_pData, other.m_pData, sizeof(double)*m_nNumColumns*m_nNumRows);
  118. }
  120. //////////////////////////////////////////////////////////////////////
  121. // 析构函数
  122. //////////////////////////////////////////////////////////////////////
  123. CMatrix::~CMatrix()
  124. {
  125. if (m_pData)
  126. {
  127. delete[] m_pData;
  128. m_pData = NULL;
  129. }
  130. }
  132. //////////////////////////////////////////////////////////////////////
  133. // 初始化函数
  134. //
  135. // 参数：
  136. // 1. int nRows - 指定的矩阵行数
  137. // 2. int nCols - 指定的矩阵列数
  138. //
  139. // 返回值：BOOL 型，初始化是否成功
  140. //////////////////////////////////////////////////////////////////////
  141. BOOL CMatrix::Init(int nRows, int nCols)
  142. {
  143. if (m_pData)
  144. {
  145. delete[] m_pData;
  146. m_pData = NULL;
  147. }
  149. m_nNumRows = nRows;
  150. m_nNumColumns = nCols;
  151. int nSize = nCols*nRows;
  152. if (nSize < 0)
  153. return FALSE;
  155. // 分配内存
  156. m_pData = new double[nSize];
  158. if (m_pData == NULL)
  159. return FALSE; // 内存分配失败
  160. if (IsBadReadPtr(m_pData, sizeof(double) * nSize))
  161. return FALSE;
  163. // 将各元素值置0
  164. memset(m_pData, 0, sizeof(double) * nSize);
  166. return TRUE;
  167. }
  169. //////////////////////////////////////////////////////////////////////
  170. // 将方阵初始化为单位矩阵
  171. //
  172. // 参数：
  173. // 1. int nSize - 方阵行列数
  174. //
  175. // 返回值：BOOL 型，初始化是否成功
  176. //////////////////////////////////////////////////////////////////////
  177. BOOL CMatrix::MakeUnitMatrix(int nSize)
  178. {
  179. if (! Init(nSize, nSize))
  180. return FALSE;
  182. for (int i=0; i<nSize; ++i)
  183. for (int j=0; j<nSize; ++j)
  184. if (i == j)
  185. SetElement(i, j, 1);
  187. return TRUE;
  188. }
  190. //////////////////////////////////////////////////////////////////////
  191. // 将字符串转化为矩阵的值
  192. //
  193. // 参数：
  194. // 1. CString s - 数字和分隔符构成的字符串
  195. // 2. const CString& sDelim - 数字之间的分隔符，默认为空格
  196. // 3. BOOL bLineBreak - 行与行之间是否有回车换行符，默认为真(有换行符)
  197. //         当该参数为FALSE时，所有元素值都在一行中输入，字符串的第一个
  198. //         数值应为矩阵的行数，第二个数值应为矩阵的列数
  199. //
  200. // 返回值：BOOL 型，转换是否成功
  201. //////////////////////////////////////////////////////////////////////
  202. BOOL CMatrix::FromString(CString s, const CString& sDelim /*= " "*/, BOOL bLineBreak /*= TRUE*/)
  203. {
  204. if (s.IsEmpty())
  205. return FALSE;
  207. // 分行处理
  208. if (bLineBreak)
  209. {
  210. CTokenizer tk(s, "rn");
  212. CStringList ListRow;
  213. CString sRow;
  214. while (tk.Next(sRow))
  215. {
  216. sRow.TrimLeft();
  217. sRow.TrimRight();
  218. if (sRow.IsEmpty())
  219. break;
  221. ListRow.AddTail(sRow);
  222. }
  224. // 行数
  225. m_nNumRows = ListRow.GetCount();
  227. sRow = ListRow.GetHead();
  228. CTokenizer tkRow(sRow, sDelim);
  229. CString sElement;
  230. // 列数
  231. m_nNumColumns = 0;
  232. while (tkRow.Next(sElement))
  233. {
  234. m_nNumColumns++;
  235. }
  237. // 初始化矩阵
  238. if (! Init(m_nNumRows, m_nNumColumns))
  239. return FALSE;
  241. // 设置值
  242. POSITION pos = ListRow.GetHeadPosition();
  243. for (int i=0; i<m_nNumRows; i++)
  244. {
  245. sRow = ListRow.GetNext(pos);
  246. int j = 0;
  247. CTokenizer tkRow(sRow, sDelim);
  248. while (tkRow.Next(sElement))
  249. {
  250. sElement.TrimLeft();
  251. sElement.TrimRight();
  252. double v = atof(sElement);
  253. SetElement(i, j++, v);
  254. }
  255. }
  257. return TRUE;
  258. }
  260. // 不分行（单行）处理
  262. CTokenizer tk(s, sDelim);
  264. CString sElement;
  266. // 行数
  267. tk.Next(sElement);
  268. sElement.TrimLeft();
  269. sElement.TrimRight();
  270. m_nNumRows = atoi(sElement);
  272. // 列数
  273. tk.Next(sElement);
  274. sElement.TrimLeft();
  275. sElement.TrimRight();
  276. m_nNumColumns = atoi(sElement);
  278. // 初始化矩阵
  279. if (! Init(m_nNumRows, m_nNumColumns))
  280. return FALSE;
  282. // 设置值
  283. int i = 0, j = 0;
  284. while (tk.Next(sElement))
  285. {
  286. sElement.TrimLeft();
  287. sElement.TrimRight();
  288. double v = atof(sElement);
  289. SetElement(i, j++, v);
  291. if (j == m_nNumColumns)
  292. {
  293. j = 0;
  294. i++;
  295. if (i == m_nNumRows)
  296. break;
  297. }
  298. }
  300. return TRUE;
  301. }
  303. //////////////////////////////////////////////////////////////////////
  304. // 将矩阵各元素的值转化为字符串
  305. //
  306. // 参数：
  307. // 1. const CString& sDelim - 数字之间的分隔符，默认为空格
  308. // 2 BOOL bLineBreak - 行与行之间是否有回车换行符，默认为真(有换行符)
  309. //
  310. // 返回值：CString 型，转换得到的字符串
  311. //////////////////////////////////////////////////////////////////////
  312. CString CMatrix::ToString(const CString& sDelim /*= " "*/, BOOL bLineBreak /*= TRUE*/) const
  313. {
  314. CString s="";
  316. for (int i=0; i<m_nNumRows; ++i)
  317. {
  318. for (int j=0; j<m_nNumColumns; ++j)
  319. {
  320. CString ss;
  321. ss.Format("%f", GetElement(i, j));
  322. s += ss;
  324. if (bLineBreak)
  325. {
  326. if (j != m_nNumColumns-1)
  327. s += sDelim;
  328. }
  329. else
  330. {
  331. if (i != m_nNumRows-1 || j != m_nNumColumns-1)
  332. s += sDelim;
  333. }
  334. }
  335. if (bLineBreak)
  336. if (i != m_nNumRows-1)
  337. s += "rn";
  338. }
  340. return s;
  341. }
  343. //////////////////////////////////////////////////////////////////////
  344. // 将矩阵指定行中各元素的值转化为字符串
  345. //
  346. // 参数：
  347. // 1. int nRow - 指定的矩阵行，nRow = 0表示第一行
  348. // 2. const CString& sDelim - 数字之间的分隔符，默认为空格
  349. //
  350. // 返回值：CString 型，转换得到的字符串
  351. //////////////////////////////////////////////////////////////////////
  352. CString CMatrix::RowToString(int nRow, const CString& sDelim /*= " "*/) const
  353. {
  354. CString s = "";
  356. if (nRow >= m_nNumRows)
  357. return s;
  359. for (int j=0; j<m_nNumColumns; ++j)
  360. {
  361. CString ss;
  362. ss.Format("%f", GetElement(nRow, j));
  363. s += ss;
  364. if (j != m_nNumColumns-1)
  365. s += sDelim;
  366. }
  368. return s;
  369. }
  371. //////////////////////////////////////////////////////////////////////
  372. // 将矩阵指定列中各元素的值转化为字符串
  373. //
  374. // 参数：
  375. // 1. int nCol - 指定的矩阵行，nCol = 0表示第一列
  376. // 2. const CString& sDelim - 数字之间的分隔符，默认为空格
  377. //
  378. // 返回值：CString 型，转换得到的字符串
  379. //////////////////////////////////////////////////////////////////////
  380. CString CMatrix::ColToString(int nCol, const CString& sDelim /*= " "*/) const
  381. {
  382. CString s = "";
  384. if (nCol >= m_nNumColumns)
  385. return s;
  387. for (int i=0; i<m_nNumRows; ++i)
  388. {
  389. CString ss;
  390. ss.Format("%f", GetElement(i, nCol));
  391. s += ss;
  392. if (i != m_nNumRows-1)
  393. s += sDelim;
  394. }
  396. return s;
  397. }
  399. //////////////////////////////////////////////////////////////////////
  400. // 设置矩阵各元素的值
  401. //
  402. // 参数：
  403. // 1. double value[] - 一维数组，长度为m_nNumColumns*m_nNumRows，存储
  404. //                     矩阵各元素的值
  405. //
  406. // 返回值：无
  407. //////////////////////////////////////////////////////////////////////
  408. void CMatrix::SetData(double value[])
  409. {
  410. // empty the memory
  411. memset(m_pData, 0, sizeof(double) * m_nNumColumns*m_nNumRows);
  412. // copy data
  413. memcpy(m_pData, value, sizeof(double)*m_nNumColumns*m_nNumRows);
  414. }
  416. //////////////////////////////////////////////////////////////////////
  417. // 设置指定元素的值
  418. //
  419. // 参数：
  420. // 1. int nRows - 指定的矩阵行数
  421. // 2. int nCols - 指定的矩阵列数
  422. // 3. double value - 指定元素的值
  423. //
  424. // 返回值：BOOL 型，说明设置是否成功
  425. //////////////////////////////////////////////////////////////////////
  426. BOOL CMatrix::SetElement(int nRow, int nCol, double value)
  427. {
  428. if (nCol < 0 || nCol >= m_nNumColumns || nRow < 0 || nRow >= m_nNumRows)
  429. return FALSE; // array bounds error
  430. if (m_pData == NULL)
  431. return FALSE; // bad pointer error
  433. m_pData[nCol + nRow * m_nNumColumns] = value;
  435. return TRUE;
  436. }
  438. //////////////////////////////////////////////////////////////////////
  439. // 设置指定元素的值
  440. //
  441. // 参数：
  442. // 1. int nRows - 指定的矩阵行数
  443. // 2. int nCols - 指定的矩阵列数
  444. //
  445. // 返回值：double 型，指定元素的值
  446. //////////////////////////////////////////////////////////////////////
  447. double CMatrix::GetElement(int nRow, int nCol) const
  448. {
  449. ASSERT(nCol >= 0 && nCol < m_nNumColumns && nRow >= 0 && nRow < m_nNumRows); // array bounds error
  450. ASSERT(m_pData); // bad pointer error
  451. return m_pData[nCol + nRow * m_nNumColumns] ;
  452. }
  454. //////////////////////////////////////////////////////////////////////
  455. // 获取矩阵的列数
  456. //
  457. // 参数：无
  458. //
  459. // 返回值：int 型，矩阵的列数
  460. //////////////////////////////////////////////////////////////////////
  461. int CMatrix::GetNumColumns() const
  462. {
  463. return m_nNumColumns;
  464. }
  466. //////////////////////////////////////////////////////////////////////
  467. // 获取矩阵的行数
  468. //
  469. // 参数：无
  470. //
  471. // 返回值：int 型，矩阵的行数
  472. //////////////////////////////////////////////////////////////////////
  473. int CMatrix::GetNumRows() const
  474. {
  475. return m_nNumRows;
  476. }
  478. //////////////////////////////////////////////////////////////////////
  479. // 获取矩阵的数据
  480. //
  481. // 参数：无
  482. //
  483. // 返回值：double型指针，指向矩阵各元素的数据缓冲区
  484. //////////////////////////////////////////////////////////////////////
  485. double* CMatrix::GetData() const
  486. {
  487. return m_pData;
  488. }
  490. //////////////////////////////////////////////////////////////////////
  491. // 获取指定行的向量
  492. //
  493. // 参数：
  494. // 1. int nRows - 指定的矩阵行数
  495. // 2.  double* pVector - 指向向量中各元素的缓冲区
  496. //
  497. // 返回值：int 型，向量中元素的个数，即矩阵的列数
  498. //////////////////////////////////////////////////////////////////////
  499. int CMatrix::GetRowVector(int nRow, double* pVector) const
  500. {
  501. if (pVector == NULL)
  502. delete pVector;
  504. pVector = new double[m_nNumColumns];
  505. ASSERT(pVector != NULL);
  507. for (int j=0; j<m_nNumColumns; ++j)
  508. pVector[j] = GetElement(nRow, j);
  510. return m_nNumColumns;
  511. }
  513. //////////////////////////////////////////////////////////////////////
  514. // 获取指定列的向量
  515. //
  516. // 参数：
  517. // 1. int nCols - 指定的矩阵列数
  518. // 2.  double* pVector - 指向向量中各元素的缓冲区
  519. //
  520. // 返回值：int 型，向量中元素的个数，即矩阵的行数
  521. //////////////////////////////////////////////////////////////////////
  522. int CMatrix::GetColVector(int nCol, double* pVector) const
  523. {
  524. if (pVector == NULL)
  525. delete pVector;
  527. pVector = new double[m_nNumRows];
  528. ASSERT(pVector != NULL);
  530. for (int i=0; i<m_nNumRows; ++i)
  531. pVector[i] = GetElement(i, nCol);
  533. return m_nNumRows;
  534. }
  536. //////////////////////////////////////////////////////////////////////
  537. // 重载运算符=，给矩阵赋值
  538. //
  539. // 参数：
  540. // 1. const CMatrix& other - 用于给矩阵赋值的源矩阵
  541. //
  542. // 返回值：CMatrix型的引用，所引用的矩阵与other相等
  543. //////////////////////////////////////////////////////////////////////
  544. CMatrix& CMatrix::operator=(const CMatrix& other)
  545. {
  546. if (&other != this)
  547. {
  548. BOOL bSuccess = Init(other.GetNumRows(), other.GetNumColumns());
  549. ASSERT(bSuccess);
  551. // copy the pointer
  552. memcpy(m_pData, other.m_pData, sizeof(double)*m_nNumColumns*m_nNumRows);
  553. }
  555. // finally return a reference to ourselves
  556. return *this ;
  557. }
  559. //////////////////////////////////////////////////////////////////////
  560. // 重载运算符==，判断矩阵是否相等
  561. //
  562. // 参数：
  563. // 1. const CMatrix& other - 用于比较的矩阵
  564. //
  565. // 返回值：BOOL 型，两个矩阵相等则为TRUE，否则为FALSE
  566. //////////////////////////////////////////////////////////////////////
  567. BOOL CMatrix::operator==(const CMatrix& other) const
  568. {
  569. // 首先检查行列数是否相等
  570. if (m_nNumColumns != other.GetNumColumns() || m_nNumRows != other.GetNumRows())
  571. return FALSE;
  573. for (int i=0; i<m_nNumRows; ++i)
  574. {
  575. for (int j=0; j<m_nNumColumns; ++j)
  576. {
  577. if (GetElement(i, j) != other.GetElement(i, j))
  578. return FALSE;
  579. }
  580. }
  582. return TRUE;
  583. }
  585. //////////////////////////////////////////////////////////////////////
  586. // 重载运算符!=，判断矩阵是否不相等
  587. //
  588. // 参数：
  589. // 1. const CMatrix& other - 用于比较的矩阵
  590. //
  591. // 返回值：BOOL 型，两个不矩阵相等则为TRUE，否则为FALSE
  592. //////////////////////////////////////////////////////////////////////
  593. BOOL CMatrix::operator!=(const CMatrix& other) const
  594. {
  595. return !(*this == other);
  596. }
  598. //////////////////////////////////////////////////////////////////////
  599. // 重载运算符+，实现矩阵的加法
  600. //
  601. // 参数：
  602. // 1. const CMatrix& other - 与指定矩阵相加的矩阵
  603. //
  604. // 返回值：CMatrix型，指定矩阵与other相加之和
  605. //////////////////////////////////////////////////////////////////////
  606. CMatrix CMatrix::operator+(const CMatrix& other) const
  607. {
  608. // 首先检查行列数是否相等
  609. ASSERT (m_nNumColumns == other.GetNumColumns() && m_nNumRows == other.GetNumRows());
  611. // 构造结果矩阵
  612. CMatrix result(*this) ; // 拷贝构造
  613. // 矩阵加法
  614. for (int i = 0 ; i < m_nNumRows ; ++i)
  615. {
  616. for (int j = 0 ; j <  m_nNumColumns; ++j)
  617. result.SetElement(i, j, result.GetElement(i, j) + other.GetElement(i, j)) ;
  618. }
  620. return result ;
  621. }
  623. //////////////////////////////////////////////////////////////////////
  624. // 重载运算符-，实现矩阵的减法
  625. //
  626. // 参数：
  627. // 1. const CMatrix& other - 与指定矩阵相减的矩阵
  628. //
  629. // 返回值：CMatrix型，指定矩阵与other相减之差
  630. //////////////////////////////////////////////////////////////////////
  631. CMatrix CMatrix::operator-(const CMatrix& other) const
  632. {
  633. // 首先检查行列数是否相等
  634. ASSERT (m_nNumColumns == other.GetNumColumns() && m_nNumRows == other.GetNumRows());
  636. // 构造目标矩阵
  637. CMatrix result(*this) ; // copy ourselves
  638. // 进行减法操作
  639. for (int i = 0 ; i < m_nNumRows ; ++i)
  640. {
  641. for (int j = 0 ; j <  m_nNumColumns; ++j)
  642. result.SetElement(i, j, result.GetElement(i, j) - other.GetElement(i, j)) ;
  643. }
  645. return result ;
  646. }
  648. //////////////////////////////////////////////////////////////////////
  649. // 重载运算符*，实现矩阵的数乘
  650. //
  651. // 参数：
  652. // 1. double value - 与指定矩阵相乘的实数
  653. //
  654. // 返回值：CMatrix型，指定矩阵与value相乘之积
  655. //////////////////////////////////////////////////////////////////////
  656. CMatrix CMatrix::operator*(double value) const
  657. {
  658. // 构造目标矩阵
  659. CMatrix result(*this) ; // copy ourselves
  660. // 进行数乘
  661. for (int i = 0 ; i < m_nNumRows ; ++i)
  662. {
  663. for (int j = 0 ; j <  m_nNumColumns; ++j)
  664. result.SetElement(i, j, result.GetElement(i, j) * value) ;
  665. }
  667. return result ;
  668. }
  670. //////////////////////////////////////////////////////////////////////
  671. // 重载运算符*，实现矩阵的乘法
  672. //
  673. // 参数：
  674. // 1. const CMatrix& other - 与指定矩阵相乘的矩阵
  675. //
  676. // 返回值：CMatrix型，指定矩阵与other相乘之积
  677. //////////////////////////////////////////////////////////////////////
  678. CMatrix CMatrix::operator*(const CMatrix& other) const
  679. {
  680. // 首先检查行列数是否符合要求
  681. ASSERT (m_nNumColumns == other.GetNumRows());
  683. // construct the object we are going to return
  684. CMatrix result(m_nNumRows, other.GetNumColumns()) ;
  686. // 矩阵乘法，即
  687. //
  688. // [A][B][C]   [G][H]     [A*G + B*I + C*K][A*H + B*J + C*L]
  689. // [D][E][F] * [I][J] =   [D*G + E*I + F*K][D*H + E*J + F*L]
  690. //             [K][L]
  691. //
  692. double value ;
  693. for (int i = 0 ; i < result.GetNumRows() ; ++i)
  694. {
  695. for (int j = 0 ; j < other.GetNumColumns() ; ++j)
  696. {
  697. value = 0.0 ;
  698. for (int k = 0 ; k < m_nNumColumns ; ++k)
  699. {
  700. value += GetElement(i, k) * other.GetElement(k, j) ;
  701. }
  703. result.SetElement(i, j, value) ;
  704. }
  705. }
  707. return result ;
  708. }
  710. //////////////////////////////////////////////////////////////////////
  711. // 复矩阵的乘法
  712. //
  713. // 参数：
  714. // 1. const CMatrix& AR - 左边复矩阵的实部矩阵
  715. // 2. const CMatrix& AI - 左边复矩阵的虚部矩阵
  716. // 3. const CMatrix& BR - 右边复矩阵的实部矩阵
  717. // 4. const CMatrix& BI - 右边复矩阵的虚部矩阵
  718. // 5. CMatrix& CR - 乘积复矩阵的实部矩阵
  719. // 6. CMatrix& CI - 乘积复矩阵的虚部矩阵
  720. //
  721. // 返回值：BOOL型，复矩阵乘法是否成功
  722. //////////////////////////////////////////////////////////////////////
  723. BOOL CMatrix::CMul(const CMatrix& AR, const CMatrix& AI, const CMatrix& BR, const CMatrix& BI, CMatrix& CR, CMatrix& CI) const
  724. {
  725. // 首先检查行列数是否符合要求
  726. if (AR.GetNumColumns() != AI.GetNumColumns() ||
  727. AR.GetNumRows() != AI.GetNumRows() ||
  728. BR.GetNumColumns() != BI.GetNumColumns() ||
  729. BR.GetNumRows() != BI.GetNumRows() ||
  730. AR.GetNumColumns() != BR.GetNumRows())
  731. return FALSE;
  733. // 构造乘积矩阵实部矩阵和虚部矩阵
  734. CMatrix mtxCR(AR.GetNumRows(), BR.GetNumColumns()), mtxCI(AR.GetNumRows(), BR.GetNumColumns());
  735. // 复矩阵相乘
  736.     for (int i=0; i<AR.GetNumRows(); ++i)
  737. {
  738.     for (int j=0; j<BR.GetNumColumns(); ++j)
  739. {
  740. double vr = 0;
  741. double vi = 0;
  742.             for (int k =0; k<AR.GetNumColumns(); ++k)
  743. {
  744.                 double p = AR.GetElement(i, k) * BR.GetElement(k, j);
  745.                 double q = AI.GetElement(i, k) * BI.GetElement(k, j);
  746.                 double s = (AR.GetElement(i, k) + AI.GetElement(i, k)) * (BR.GetElement(k, j) + BI.GetElement(k, j));
  747.                 vr += p - q;
  748.                 vi += s - p - q;
  749. }
  750.             mtxCR.SetElement(i, j, vr);
  751.             mtxCI.SetElement(i, j, vi);
  752.         }
  753. }
  755. CR = mtxCR;
  756. CI = mtxCI;
  758. return TRUE;
  759. }
  761. //////////////////////////////////////////////////////////////////////
  762. // 矩阵的转置
  763. //
  764. // 参数：无
  765. //
  766. // 返回值：CMatrix型，指定矩阵转置矩阵
  767. //////////////////////////////////////////////////////////////////////
  768. CMatrix CMatrix::Transpose() const
  769. {
  770. // 构造目标矩阵
  771. CMatrix Trans(m_nNumColumns, m_nNumRows);
  773. // 转置各元素
  774. for (int i = 0 ; i < m_nNumRows ; ++i)
  775. {
  776. for (int j = 0 ; j < m_nNumColumns ; ++j)
  777. Trans.SetElement(j, i, GetElement(i, j)) ;
  778. }
  780. return Trans;
  781. }
  783. //////////////////////////////////////////////////////////////////////
  784. // 实矩阵求逆的全选主元高斯－约当法
  785. //
  786. // 参数：无
  787. //
  788. // 返回值：BOOL型，求逆是否成功
  789. //////////////////////////////////////////////////////////////////////
  790. BOOL CMatrix::InvertGaussJordan()
  791. {
  792. int *pnRow, *pnCol,i,j,k,l,u,v;
  793.     double d = 0, p = 0;
  795. // 分配内存
  796.     pnRow = new int[m_nNumColumns];
  797.     pnCol = new int[m_nNumColumns];
  798. if (pnRow == NULL || pnCol == NULL)
  799. return FALSE;
  801. // 消元
  802.     for (k=0; k<=m_nNumColumns-1; k++)
  803.     { 
  804. d=0.0;
  805.         for (i=k; i<=m_nNumColumns-1; i++)
  806. {
  807. for (j=k; j<=m_nNumColumns-1; j++)
  808. { 
  809. l=i*m_nNumColumns+j; p=fabs(m_pData[l]);
  810. if (p>d) 
  811. { 
  812. d=p; 
  813. pnRow[k]=i; 
  814. pnCol[k]=j;
  815. }
  816. }
  817. }
  818.         
  819. // 失败
  820. if (d == 0.0)
  821. {
  822. delete[] pnRow;
  823. delete[] pnCol;
  824. return FALSE;
  825. }
  827.         if (pnRow[k] != k)
  828. {
  829. for (j=0; j<=m_nNumColumns-1; j++)
  830. { 
  831. u=k*m_nNumColumns+j; 
  832. v=pnRow[k]*m_nNumColumns+j;
  833. p=m_pData[u]; 
  834. m_pData[u]=m_pData[v]; 
  835. m_pData[v]=p;
  836. }
  837. }
  838.         
  839. if (pnCol[k] != k)
  840. {
  841. for (i=0; i<=m_nNumColumns-1; i++)
  842.             { 
  843. u=i*m_nNumColumns+k; 
  844. v=i*m_nNumColumns+pnCol[k];
  845. p=m_pData[u]; 
  846. m_pData[u]=m_pData[v]; 
  847. m_pData[v]=p;
  848.             }
  849. }
  851.         l=k*m_nNumColumns+k;
  852.         m_pData[l]=1.0/m_pData[l];
  853.         for (j=0; j<=m_nNumColumns-1; j++)
  854. {
  855. if (j != k)
  856.             { 
  857. u=k*m_nNumColumns+j; 
  858. m_pData[u]=m_pData[u]*m_pData[l];
  859. }
  860. }
  862.         for (i=0; i<=m_nNumColumns-1; i++)
  863. {
  864. if (i!=k)
  865. {
  866. for (j=0; j<=m_nNumColumns-1; j++)
  867. {
  868. if (j!=k)
  869. { 
  870. u=i*m_nNumColumns+j;
  871. m_pData[u]=m_pData[u]-m_pData[i*m_nNumColumns+k]*m_pData[k*m_nNumColumns+j];
  872. }
  873.                 }
  874. }
  875. }
  877.         for (i=0; i<=m_nNumColumns-1; i++)
  878. {
  879. if (i!=k)
  880.             { 
  881. u=i*m_nNumColumns+k; 
  882. m_pData[u]=-m_pData[u]*m_pData[l];
  883. }
  884. }
  885.     }
  887.     // 调整恢复行列次序
  888.     for (k=m_nNumColumns-1; k>=0; k--)
  889.     { 
  890. if (pnCol[k]!=k)
  891. {
  892. for (j=0; j<=m_nNumColumns-1; j++)
  893.             { 
  894. u=k*m_nNumColumns+j; 
  895. v=pnCol[k]*m_nNumColumns+j;
  896. p=m_pData[u]; 
  897. m_pData[u]=m_pData[v]; 
  898. m_pData[v]=p;
  899.             }
  900. }
  902.         if (pnRow[k]!=k)
  903. {
  904. for (i=0; i<=m_nNumColumns-1; i++)
  905.             { 
  906. u=i*m_nNumColumns+k; 
  907. v=i*m_nNumColumns+pnRow[k];
  908. p=m_pData[u]; 
  909. m_pData[u]=m_pData[v]; 
  910. m_pData[v]=p;
  911.             }
  912. }
  913.     }
  915. // 清理内存
  916. delete[] pnRow;
  917. delete[] pnCol;
  919. // 成功返回
  920. return TRUE;
  921. }
  923. //////////////////////////////////////////////////////////////////////
  924. // 复矩阵求逆的全选主元高斯－约当法
  925. //
  926. // 参数：
  927. // 1. CMatrix& mtxImag - 复矩阵的虚部矩阵，当前矩阵为复矩阵的实部
  928. //
  929. // 返回值：BOOL型，求逆是否成功
  930. //////////////////////////////////////////////////////////////////////
  931. BOOL CMatrix::InvertGaussJordan(CMatrix& mtxImag)
  932. {
  933. int *pnRow,*pnCol,i,j,k,l,u,v,w;
  934.     double p,q,s,t,d,b;
  936. // 分配内存
  937.     pnRow = new int[m_nNumColumns];
  938.     pnCol = new int[m_nNumColumns];
  939. if (pnRow == NULL || pnCol == NULL)
  940. return FALSE;
  942. // 消元
  943.     for (k=0; k<=m_nNumColumns-1; k++)
  944.     { 
  945. d=0.0;
  946.         for (i=k; i<=m_nNumColumns-1; i++)
  947. {
  948. for (j=k; j<=m_nNumColumns-1; j++)
  949. { 
  950. u=i*m_nNumColumns+j;
  951. p=m_pData[u]*m_pData[u]+mtxImag.m_pData[u]*mtxImag.m_pData[u];
  952. if (p>d) 
  953. { 
  954. d=p; 
  955. pnRow[k]=i; 
  956. pnCol[k]=j;
  957. }
  958. }
  959. }
  961. // 失败
  962.         if (d == 0.0)
  963.         { 
  964. delete[] pnRow;
  965. delete[] pnCol;
  966.             return(0);
  967.         }
  969.         if (pnRow[k]!=k)
  970. {
  971. for (j=0; j<=m_nNumColumns-1; j++)
  972.             { 
  973. u=k*m_nNumColumns+j; 
  974. v=pnRow[k]*m_nNumColumns+j;
  975. t=m_pData[u]; 
  976. m_pData[u]=m_pData[v]; 
  977. m_pData[v]=t;
  978. t=mtxImag.m_pData[u]; 
  979. mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
  980. mtxImag.m_pData[v]=t;
  981.             }
  982. }
  984.         if (pnCol[k]!=k)
  985. {
  986. for (i=0; i<=m_nNumColumns-1; i++)
  987.             { 
  988. u=i*m_nNumColumns+k; 
  989. v=i*m_nNumColumns+pnCol[k];
  990. t=m_pData[u]; 
  991. m_pData[u]=m_pData[v]; 
  992. m_pData[v]=t;
  993. t=mtxImag.m_pData[u]; 
  994. mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
  995. mtxImag.m_pData[v]=t;
  996.             }
  997. }
  999.         l=k*m_nNumColumns+k;
  1000.         m_pData[l]=m_pData[l]/d; mtxImag.m_pData[l]=-mtxImag.m_pData[l]/d;
  1001.         for (j=0; j<=m_nNumColumns-1; j++)
  1002. {
  1003. if (j!=k)
  1004.             { 
  1005. u=k*m_nNumColumns+j;
  1006. p=m_pData[u]*m_pData[l]; 
  1007. q=mtxImag.m_pData[u]*mtxImag.m_pData[l];
  1008. s=(m_pData[u]+mtxImag.m_pData[u])*(m_pData[l]+mtxImag.m_pData[l]);
  1009. m_pData[u]=p-q; 
  1010. mtxImag.m_pData[u]=s-p-q;
  1011.             }
  1012. }
  1014.         for (i=0; i<=m_nNumColumns-1; i++)
  1015. {
  1016. if (i!=k)
  1017.             { 
  1018. v=i*m_nNumColumns+k;
  1019. for (j=0; j<=m_nNumColumns-1; j++)
  1020. {
  1021. if (j!=k)
  1022. { 
  1023. u=k*m_nNumColumns+j;  
  1024. w=i*m_nNumColumns+j;
  1025. p=m_pData[u]*m_pData[v]; 
  1026. q=mtxImag.m_pData[u]*mtxImag.m_pData[v];
  1027. s=(m_pData[u]+mtxImag.m_pData[u])*(m_pData[v]+mtxImag.m_pData[v]);
  1028. t=p-q; 
  1029. b=s-p-q;
  1030. m_pData[w]=m_pData[w]-t;
  1031. mtxImag.m_pData[w]=mtxImag.m_pData[w]-b;
  1032. }
  1033. }
  1034.             }
  1035. }
  1037.         for (i=0; i<=m_nNumColumns-1; i++)
  1038. {
  1039. if (i!=k)
  1040.             { 
  1041. u=i*m_nNumColumns+k;
  1042. p=m_pData[u]*m_pData[l]; 
  1043. q=mtxImag.m_pData[u]*mtxImag.m_pData[l];
  1044. s=(m_pData[u]+mtxImag.m_pData[u])*(m_pData[l]+mtxImag.m_pData[l]);
  1045. m_pData[u]=q-p; 
  1046. mtxImag.m_pData[u]=p+q-s;
  1047.             }
  1048. }
  1049.     }
  1051.     // 调整恢复行列次序
  1052.     for (k=m_nNumColumns-1; k>=0; k--)
  1053.     { 
  1054. if (pnCol[k]!=k)
  1055. {
  1056. for (j=0; j<=m_nNumColumns-1; j++)
  1057.             { 
  1058. u=k*m_nNumColumns+j; 
  1059. v=pnCol[k]*m_nNumColumns+j;
  1060. t=m_pData[u]; 
  1061. m_pData[u]=m_pData[v]; 
  1062. m_pData[v]=t;
  1063. t=mtxImag.m_pData[u]; 
  1064. mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
  1065. mtxImag.m_pData[v]=t;
  1066.             }
  1067. }
  1069.         if (pnRow[k]!=k)
  1070. {
  1071. for (i=0; i<=m_nNumColumns-1; i++)
  1072.             { 
  1073. u=i*m_nNumColumns+k; 
  1074. v=i*m_nNumColumns+pnRow[k];
  1075. t=m_pData[u]; 
  1076. m_pData[u]=m_pData[v]; 
  1077. m_pData[v]=t;
  1078. t=mtxImag.m_pData[u]; 
  1079. mtxImag.m_pData[u]=mtxImag.m_pData[v]; 
  1080. mtxImag.m_pData[v]=t;
  1081.             }
  1082. }
  1083.     }
  1085. // 清理内存
  1086. delete[] pnRow;
  1087. delete[] pnCol;
  1089. // 成功返回
  1090. return TRUE;
  1091. }
  1093. //////////////////////////////////////////////////////////////////////
  1094. // 对称正定矩阵的求逆
  1095. //
  1096. // 参数：无
  1097. //
  1098. // 返回值：BOOL型，求逆是否成功
  1099. //////////////////////////////////////////////////////////////////////
  1100. BOOL CMatrix::InvertSsgj()
  1101. { 
  1102. int i, j ,k, m;
  1103.     double w, g, *pTmp;
  1105. // 临时内存
  1106.     pTmp = new double[m_nNumColumns];
  1108. // 逐列处理
  1109.     for (k=0; k<=m_nNumColumns-1; k++)
  1110.     { 
  1111. w=m_pData[0];
  1112.         if (w == 0.0)
  1113.         { 
  1114. delete[] pTmp;
  1115. return FALSE;
  1116. }
  1118.         m=m_nNumColumns-k-1;
  1119.         for (i=1; i<=m_nNumColumns-1; i++)
  1120.         { 
  1121. g=m_pData[i*m_nNumColumns]; 
  1122. pTmp[i]=g/w;
  1123.             if (i<=m) 
  1124. pTmp[i]=-pTmp[i];
  1125.             for (j=1; j<=i; j++)
  1126.               m_pData[(i-1)*m_nNumColumns+j-1]=m_pData[i*m_nNumColumns+j]+g*pTmp[j];
  1127.         }
  1129.         m_pData[m_nNumColumns*m_nNumColumns-1]=1.0/w;
  1130.         for (i=1; i<=m_nNumColumns-1; i++)
  1131. m_pData[(m_nNumColumns-1)*m_nNumColumns+i-1]=pTmp[i];
  1132.     }
  1134. // 行列调整
  1135.     for (i=0; i<=m_nNumColumns-2; i++)
  1136. for (j=i+1; j<=m_nNumColumns-1; j++)
  1137. m_pData[i*m_nNumColumns+j]=m_pData[j*m_nNumColumns+i];
  1139. // 临时内存清理
  1140. delete[] pTmp;
  1142. return TRUE;
  1143. }
  1145. //////////////////////////////////////////////////////////////////////
  1146. // 托伯利兹矩阵求逆的埃兰特方法
  1147. //
  1148. // 参数：无
  1149. //
  1150. // 返回值：BOOL型，求逆是否成功
  1151. //////////////////////////////////////////////////////////////////////
  1152. BOOL CMatrix::InvertTrench()
  1153. { 
  1154. int i,j,k;
  1155.     double a,s,*t,*tt,*c,*r,*p;
  1157. // 上三角元素
  1158. t = new double[m_nNumColumns];
  1159. // 下三角元素
  1160. tt = new double[m_nNumColumns];
  1162. // 上、下三角元素赋值
  1163. for (i=0; i<m_nNumColumns; ++i)
  1164. {
  1165. t[i] = GetElement(0, i);
  1166.     tt[i] = GetElement(i, 0);
  1167. }
  1169. // 临时缓冲区
  1170. c = new double[m_nNumColumns];
  1171. r = new double[m_nNumColumns];
  1172. p = new double[m_nNumColumns];
  1174. // 非Toeplitz矩阵，返回
  1175.     if (t[0] == 0.0)
  1176.     { 
  1177. delete[] t;
  1178. delete[] tt;
  1179. delete[] c;
  1180. delete[] r;
  1181. delete[] p;
  1182. return FALSE;
  1183.     }
  1185.     a=t[0]; 
  1186. c[0]=tt[1]/t[0]; 
  1187. r[0]=t[1]/t[0];
  1189.     for (k=0; k<=m_nNumColumns-3; k++)
  1190.     { 
  1191. s=0.0;
  1192.         for (j=1; j<=k+1; j++)
  1193. s=s+c[k+1-j]*tt[j];
  1195.         s=(s-tt[k+2])/a;
  1196. for (i=0; i<=k; i++)
  1197. p[i]=c[i]+s*r[k-i];
  1199.         c[k+1]=-s;
  1200.         s=0.0;
  1201.         for (j=1; j<=k+1; j++)
  1202.           s=s+r[k+1-j]*t[j];
  1203.         
  1204. s=(s-t[k+2])/a;
  1205.         for (i=0; i<=k; i++)
  1206.         { 
  1207. r[i]=r[i]+s*c[k-i];
  1208.             c[k-i]=p[k-i];
  1209.         }
  1211.         r[k+1]=-s;
  1212. a=0.0;
  1213.         for (j=1; j<=k+2; j++)
  1214.           a=a+t[j]*c[j-1];
  1216.         a=t[0]-a;
  1218. // 求解失败
  1219.         if (a == 0.0)
  1220. { 
  1221. delete[] t;
  1222. delete[] tt;
  1223. delete[] c;
  1224. delete[] r;
  1225. delete[] p;
  1226. return FALSE;
  1227. }
  1228.     }
  1230.     m_pData[0]=1.0/a;
  1231.     for (i=0; i<=m_nNumColumns-2; i++)
  1232.     { 
  1233. k=i+1; 
  1234. j=(i+1)*m_nNumColumns;
  1235.         m_pData[k]=-r[i]/a; 
  1236. m_pData[j]=-c[i]/a;
  1237.     }
  1239.    for (i=0; i<=m_nNumColumns-2; i++)
  1240. {
  1241. for (j=0; j<=m_nNumColumns-2; j++)
  1242. { 
  1243. k=(i+1)*m_nNumColumns+j+1;
  1244. m_pData[k]=m_pData[i*m_nNumColumns+j]-c[i]*m_pData[j+1];
  1245. m_pData[k]=m_pData[k]+c[m_nNumColumns-j-2]*m_pData[m_nNumColumns-i-1];
  1246. }
  1247. }
  1249.     // 临时内存清理
  1250. delete[] t;
  1251. delete[] tt;
  1252. delete[] c;
  1253. delete[] r;
  1254. delete[] p;
  1256. return TRUE;
  1257. }
  1258.                                                
  1259. //////////////////////////////////////////////////////////////////////
  1260. // 求行列式值的全选主元高斯消去法
  1261. //
  1262. // 参数：无
  1263. //
  1264. // 返回值：double型，行列式的值
  1265. //////////////////////////////////////////////////////////////////////
  1266. double CMatrix::DetGauss()
  1267. { 
  1268. int i,j,k,is,js,l,u,v;
  1269.     double f,det,q,d;
  1270.     
  1271. // 初值
  1272. f=1.0; 
  1273. det=1.0;
  1274.     
  1275. // 消元
  1276. for (k=0; k<=m_nNumColumns-2; k++)
  1277.     { 
  1278. q=0.0;
  1279.         for (i=k; i<=m_nNumColumns-1; i++)
  1280. {
  1281. for (j=k; j<=m_nNumColumns-1; j++)
  1282. { 
  1283. l=i*m_nNumColumns+j; 
  1284. d=fabs(m_pData[l]);
  1285. if (d>q) 
  1286. { 
  1287. q=d; 
  1288. is=i; 
  1289. js=j;
  1290. }
  1291. }
  1292. }
  1294.         if (q == 0.0)
  1295.         { 
  1296. det=0.0; 
  1297. return(det);
  1298. }
  1299.         
  1300. if (is!=k)
  1301.         { 
  1302. f=-f;
  1303.             for (j=k; j<=m_nNumColumns-1; j++)
  1304.             { 
  1305. u=k*m_nNumColumns+j; 
  1306. v=is*m_nNumColumns+j;
  1307.                 d=m_pData[u]; 
  1308. m_pData[u]=m_pData[v]; 
  1309. m_pData[v]=d;
  1310.             }
  1311.         }
  1312.         
  1313. if (js!=k)
  1314.         { 
  1315. f=-f;
  1316.             for (i=k; i<=m_nNumColumns-1; i++)
  1317.             {
  1318. u=i*m_nNumColumns+js; 
  1319. v=i*m_nNumColumns+k;
  1320.                 d=m_pData[u]; 
  1321. m_pData[u]=m_pData[v]; 
  1322. m_pData[v]=d;
  1323.             }
  1324.         }
  1326.         l=k*m_nNumColumns+k;
  1327.         det=det*m_pData[l];
  1328.         for (i=k+1; i<=m_nNumColumns-1; i++)
  1329.         { 
  1330. d=m_pData[i*m_nNumColumns+k]/m_pData[l];
  1331.             for (j=k+1; j<=m_nNumColumns-1; j++)
  1332.             { 
  1333. u=i*m_nNumColumns+j;
  1334.                 m_pData[u]=m_pData[u]-d*m_pData[k*m_nNumColumns+j];
  1335.             }
  1336.         }
  1337.     }
  1338.     
  1339. // 求值
  1340. det=f*det*m_pData[m_nNumColumns*m_nNumColumns-1];
  1342.     return(det);
  1343. }
  1345. //////////////////////////////////////////////////////////////////////
  1346. // 求矩阵秩的全选主元高斯消去法
  1347. //
  1348. // 参数：无
  1349. //
  1350. // 返回值：int型，矩阵的秩
  1351. //////////////////////////////////////////////////////////////////////
  1352. int CMatrix::RankGauss()
  1353. { 
  1354. int i,j,k,nn,is,js,l,ll,u,v;
  1355.     double q,d;
  1356.     
  1357. // 秩小于等于行列数
  1358. nn = m_nNumRows;
  1359.     if (m_nNumRows >= m_nNumColumns) 
  1360. nn = m_nNumColumns;
  1362.     k=0;
  1364. // 消元求解
  1365.     for (l=0; l<=nn-1; l++)
  1366.     { 
  1367. q=0.0;
  1368.         for (i=l; i<=m_nNumRows-1; i++)
  1369. {
  1370. for (j=l; j<=m_nNumColumns-1; j++)
  1371. { 
  1372. ll=i*m_nNumColumns+j; 
  1373. d=fabs(m_pData[ll]);
  1374. if (d>q) 
  1375. { 
  1376. q=d; 
  1377. is=i; 
  1378. js=j;
  1379. }
  1380. }
  1381. }
  1383.         if (q == 0.0) 
  1384. return(k);
  1386.         k=k+1;
  1387.         if (is!=l)
  1388.         { 
  1389. for (j=l; j<=m_nNumColumns-1; j++)
  1390.             { 
  1391. u=l*m_nNumColumns+j; 
  1392. v=is*m_nNumColumns+j;
  1393.                 d=m_pData[u]; 
  1394. m_pData[u]=m_pData[v]; 
  1395. m_pData[v]=d;
  1396.             }
  1397.         }
  1398.         if (js!=l)
  1399.         { 
  1400. for (i=l; i<=m_nNumRows-1; i++)
  1401.             { 
  1402. u=i*m_nNumColumns+js; 
  1403. v=i*m_nNumColumns+l;
  1404.                 d=m_pData[u]; 
  1405. m_pData[u]=m_pData[v]; 
  1406. m_pData[v]=d;
  1407.             }
  1408.         }
  1409.         
  1410. ll=l*m_nNumColumns+l;
  1411.         for (i=l+1; i<=m_nNumColumns-1; i++)
  1412.         { 
  1413. d=m_pData[i*m_nNumColumns+l]/m_pData[ll];
  1414.             for (j=l+1; j<=m_nNumColumns-1; j++)
  1415.             { 
  1416. u=i*m_nNumColumns+j;
  1417.                 m_pData[u]=m_pData[u]-d*m_pData[l*m_nNumColumns+j];
  1418.             }
  1419.         }
  1420.     }
  1421.     
  1422. return(k);
  1423. }
  1425. //////////////////////////////////////////////////////////////////////
  1426. // 对称正定矩阵的乔里斯基分解与行列式的求值
  1427. //
  1428. // 参数：
  1429. // 1. double* dblDet - 返回行列式的值
  1430. //
  1431. // 返回值：BOOL型，求解是否成功
  1432. //////////////////////////////////////////////////////////////////////
  1433. BOOL CMatrix::DetCholesky(double* dblDet)
  1434. { 
  1435. int i,j,k,u,l;
  1436.     double d;
  1437.     
  1438. // 不满足求解要求
  1439. if (m_pData[0] <= 0.0)
  1440. return FALSE;
  1442. // 乔里斯基分解
  1444.     m_pData[0]=sqrt(m_pData[0]);
  1445.     d=m_pData[0];
  1447.     for (i=1; i<=m_nNumColumns-1; i++)
  1448.     { 
  1449. u=i*m_nNumColumns; 
  1450. m_pData[u]=m_pData[u]/m_pData[0];
  1451. }
  1452.     
  1453. for (j=1; j<=m_nNumColumns-1; j++)
  1454.     { 
  1455. l=j*m_nNumColumns+j;
  1456.         for (k=0; k<=j-1; k++)
  1457.         { 
  1458. u=j*m_nNumColumns+k; 
  1459. m_pData[l]=m_pData[l]-m_pData[u]*m_pData[u];
  1460. }
  1461.         
  1462. if (m_pData[l] <= 0.0)
  1463. return FALSE;
  1465.         m_pData[l]=sqrt(m_pData[l]);
  1466.         d=d*m_pData[l];
  1467.         
  1468. for (i=j+1; i<=m_nNumColumns-1; i++)
  1469.         { 
  1470. u=i*m_nNumColumns+j;
  1471.             for (k=0; k<=j-1; k++)
  1472. m_pData[u]=m_pData[u]-m_pData[i*m_nNumColumns+k]*m_pData[j*m_nNumColumns+k];
  1473.             
  1474. m_pData[u]=m_pData[u]/m_pData[l];
  1475.         }
  1476.     }
  1477.     
  1478. // 行列式求值
  1479. *dblDet=d*d;
  1480.     
  1481. // 下三角矩阵
  1482.     for (i=0; i<=m_nNumColumns-2; i++)
  1483. for (j=i+1; j<=m_nNumColumns-1; j++)
  1484. m_pData[i*m_nNumColumns+j]=0.0;
  1486. return TRUE;
  1487. }
  1489. //////////////////////////////////////////////////////////////////////
  1490. // 矩阵的三角分解，分解成功后，原矩阵将成为Q矩阵
  1491. //
  1492. // 参数：
  1493. // 1. CMatrix& mtxL - 返回分解后的L矩阵
  1494. // 2. CMatrix& mtxU - 返回分解后的U矩阵
  1495. //
  1496. // 返回值：BOOL型，求解是否成功
  1497. //////////////////////////////////////////////////////////////////////
  1498. BOOL CMatrix::SplitLU(CMatrix& mtxL, CMatrix& mtxU)
  1499. { 
  1500. int i,j,k,w,v,ll;
  1501.     
  1502. // 初始化结果矩阵
  1503. if (! mtxL.Init(m_nNumColumns, m_nNumColumns) ||
  1504. ! mtxU.Init(m_nNumColumns, m_nNumColumns))
  1505. return FALSE;
  1507. for (k=0; k<=m_nNumColumns-2; k++)
  1508.     { 
  1509. ll=k*m_nNumColumns+k;
  1510. if (m_pData[ll] == 0.0)
  1511. return FALSE;
  1513.         for (i=k+1; i<=m_nNumColumns-1; i++)
  1514. { 
  1515. w=i*m_nNumColumns+k; 
  1516. m_pData[w]=m_pData[w]/m_pData[ll];
  1517. }
  1519.         for (i=k+1; i<=m_nNumColumns-1; i++)
  1520.         { 
  1521. w=i*m_nNumColumns+k;
  1522.             for (j=k+1; j<=m_nNumColumns-1; j++)
  1523.             { 
  1524. v=i*m_nNumColumns+j;
  1525.                 m_pData[v]=m_pData[v]-m_pData[w]*m_pData[k*m_nNumColumns+j];
  1526.             }
  1527.         }
  1528.     }
  1529.     
  1530. for (i=0; i<=m_nNumColumns-1; i++)
  1531.     {
  1532. for (j=0; j<i; j++)
  1533.         { 
  1534. w=i*m_nNumColumns+j; 
  1535. mtxL.m_pData[w]=m_pData[w]; 
  1536. mtxU.m_pData[w]=0.0;
  1537. }
  1539.         w=i*m_nNumColumns+i;
  1540.         mtxL.m_pData[w]=1.0; 
  1541. mtxU.m_pData[w]=m_pData[w];
  1542.         
  1543. for (j=i+1; j<=m_nNumColumns-1; j++)
  1544.         { 
  1545. w=i*m_nNumColumns+j; 
  1546. mtxL.m_pData[w]=0.0; 
  1547. mtxU.m_pData[w]=m_pData[w];
  1548. }
  1549.     }
  1551. return TRUE;
  1552. }
  1554. //////////////////////////////////////////////////////////////////////
  1555. // 一般实矩阵的QR分解，分解成功后，原矩阵将成为R矩阵
  1556. //
  1557. // 参数：
  1558. // 1. CMatrix& mtxQ - 返回分解后的Q矩阵
  1559. //
  1560. // 返回值：BOOL型，求解是否成功
  1561. //////////////////////////////////////////////////////////////////////
  1562. BOOL CMatrix::SplitQR(CMatrix& mtxQ)
  1563. { 
  1564. int i,j,k,l,nn,p,jj;
  1565.     double u,alpha,w,t;
  1566.     
  1567. if (m_nNumRows < m_nNumColumns)
  1568. return FALSE;
  1570. // 初始化Q矩阵
  1571. if (! mtxQ.Init(m_nNumRows, m_nNumRows))
  1572. return FALSE;
  1574. // 对角线元素单位化
  1575.     for (i=0; i<=m_nNumRows-1; i++)
  1576. {
  1577. for (j=0; j<=m_nNumRows-1; j++)
  1578. { 
  1579. l=i*m_nNumRows+j; 
  1580. mtxQ.m_pData[l]=0.0;
  1581. if (i==j) 
  1582. mtxQ.m_pData[l]=1.0;
  1583. }
  1584. }
  1586. // 开始分解
  1588.     nn=m_nNumColumns;
  1589.     if (m_nNumRows == m_nNumColumns) 
  1590. nn=m_nNumRows-1;
  1592.     for (k=0; k<=nn-1; k++)
  1593.     { 
  1594. u=0.0; 
  1595. l=k*m_nNumColumns+k;
  1596.         for (i=k; i<=m_nNumRows-1; i++)
  1597.         { 
  1598. w=fabs(m_pData[i*m_nNumColumns+k]);
  1599.             if (w>u) 
  1600. u=w;
  1601.         }
  1602.         
  1603. alpha=0.0;
  1604.         for (i=k; i<=m_nNumRows-1; i++)
  1605.         { 
  1606. t=m_pData[i*m_nNumColumns+k]/u; 
  1607. alpha=alpha+t*t;
  1608. }
  1610.         if (m_pData[l]>0.0) 
  1611. u=-u;
  1613.         alpha=u*sqrt(alpha);
  1614.         if (alpha == 0.0)
  1615. return FALSE;
  1617.         u=sqrt(2.0*alpha*(alpha-m_pData[l]));
  1618.         if ((u+1.0)!=1.0)
  1619.         { 
  1620. m_pData[l]=(m_pData[l]-alpha)/u;
  1621.             for (i=k+1; i<=m_nNumRows-1; i++)
  1622.             { 
  1623. p=i*m_nNumColumns+k; 
  1624. m_pData[p]=m_pData[p]/u;
  1625. }
  1626.             
  1627. for (j=0; j<=m_nNumRows-1; j++)
  1628.             { 
  1629. t=0.0;
  1630.                 for (jj=k; jj<=m_nNumRows-1; jj++)
  1631. t=t+m_pData[jj*m_nNumColumns+k]*mtxQ.m_pData[jj*m_nNumRows+j];
  1633.                 for (i=k; i<=m_nNumRows-1; i++)
  1634.                 { 
  1635. p=i*m_nNumRows+j; 
  1636. mtxQ.m_pData[p]=mtxQ.m_pData[p]-2.0*t*m_pData[i*m_nNumColumns+k];
  1637. }
  1638.             }
  1639.             
  1640. for (j=k+1; j<=m_nNumColumns-1; j++)
  1641.             { 
  1642. t=0.0;
  1643.                 
  1644. for (jj=k; jj<=m_nNumRows-1; jj++)
  1645. t=t+m_pData[jj*m_nNumColumns+k]*m_pData[jj*m_nNumColumns+j];
  1646.                 
  1647. for (i=k; i<=m_nNumRows-1; i++)
  1648.                 { 
  1649. p=i*m_nNumColumns+j; 
  1650. m_pData[p]=m_pData[p]-2.0*t*m_pData[i*m_nNumColumns+k];
  1651. }
  1652.             }
  1653.             
  1654. m_pData[l]=alpha;
  1655.             for (i=k+1; i<=m_nNumRows-1; i++)
  1656. m_pData[i*m_nNumColumns+k]=0.0;
  1657.         }
  1658.     }
  1659.     
  1660. // 调整元素
  1661. for (i=0; i<=m_nNumRows-2; i++)
  1662. {
  1663. for (j=i+1; j<=m_nNumRows-1;j++)
  1664. { 
  1665. p=i*m_nNumRows+j; 
  1666. l=j*m_nNumRows+i;
  1667. t=mtxQ.m_pData[p]; 
  1668. mtxQ.m_pData[p]=mtxQ.m_pData[l]; 
  1669. mtxQ.m_pData[l]=t;
  1670. }
  1671. }
  1673. return TRUE;
  1674. }
  1676. //////////////////////////////////////////////////////////////////////
  1677. // 一般实矩阵的奇异值分解，分解成功后，原矩阵对角线元素就是矩阵的奇异值
  1678. //
  1679. // 参数：
  1680. // 1. CMatrix& mtxU - 返回分解后的U矩阵
  1681. // 2. CMatrix& mtxV - 返回分解后的V矩阵
  1682. // 3. double eps - 计算精度，默认值为0.000001
  1683. //
  1684. // 返回值：BOOL型，求解是否成功
  1685. //////////////////////////////////////////////////////////////////////
  1686. BOOL CMatrix::SplitUV(CMatrix& mtxU, CMatrix& mtxV, double eps /*= 0.000001*/)
  1687. { 
  1688. int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks;
  1689.     double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];
  1690.     double *s,*e,*w;
  1692. int m = m_nNumRows;
  1693. int n = m_nNumColumns;
  1695. // 初始化U, V矩阵
  1696. if (! mtxU.Init(m, m) || ! mtxV.Init(n, n))
  1697. return FALSE;
  1699. // 临时缓冲区
  1700. int ka = max(m, n) + 1;
  1701.     s = new double[ka];
  1702.     e = new double[ka];
  1703.     w = new double[ka];
  1705. // 指定迭代次数为60
  1706.     it=60; 
  1707. k=n;
  1709.     if (m-1<n) 
  1710. k=m-1;
  1712.     l=m;
  1713.     if (n-2<m) 
  1714. l=n-2;
  1715.     if (l<0) 
  1716. l=0;
  1718. // 循环迭代计算
  1719.     ll=k;
  1720.     if (l>k) 
  1721. ll=l;
  1722.     if (ll>=1)
  1723.     { 
  1724. for (kk=1; kk<=ll; kk++)
  1725.         { 
  1726. if (kk<=k)
  1727.             { 
  1728. d=0.0;
  1729.                 for (i=kk; i<=m; i++)
  1730.                 { 
  1731. ix=(i-1)*n+kk-1; 
  1732. d=d+m_pData[ix]*m_pData[ix];
  1733. }
  1735.                 s[kk-1]=sqrt(d);
  1736.                 if (s[kk-1]!=0.0)
  1737.                 { 
  1738. ix=(kk-1)*n+kk-1;
  1739.                     if (m_pData[ix]!=0.0)
  1740.                     { 
  1741. s[kk-1]=fabs(s[kk-1]);
  1742.                         if (m_pData[ix]<0.0) 
  1743. s[kk-1]=-s[kk-1];
  1744.                     }
  1745.                     
  1746. for (i=kk; i<=m; i++)
  1747.                     { 
  1748. iy=(i-1)*n+kk-1;
  1749.                         m_pData[iy]=m_pData[iy]/s[kk-1];
  1750.                     }
  1751.                     
  1752. m_pData[ix]=1.0+m_pData[ix];
  1753.                 }
  1754.                 
  1755. s[kk-1]=-s[kk-1];
  1756.             }
  1757.             
  1758. if (n>=kk+1)
  1759.             { 
  1760. for (j=kk+1; j<=n; j++)
  1761.                 { 
  1762. if ((kk<=k)&&(s[kk-1]!=0.0))
  1763.                     { 
  1764. d=0.0;
  1765.                         for (i=kk; i<=m; i++)
  1766.                         { 
  1767. ix=(i-1)*n+kk-1;
  1768.                             iy=(i-1)*n+j-1;
  1769.                             d=d+m_pData[ix]*m_pData[iy];
  1770.                         }
  1771.                         
  1772. d=-d/m_pData[(kk-1)*n+kk-1];
  1773.                         for (i=kk; i<=m; i++)
  1774.                         { 
  1775. ix=(i-1)*n+j-1;
  1776.                             iy=(i-1)*n+kk-1;
  1777.                             m_pData[ix]=m_pData[ix]+d*m_pData[iy];
  1778.                         }
  1779.                     }
  1780.                     
  1781. e[j-1]=m_pData[(kk-1)*n+j-1];
  1782.                 }
  1783.             }
  1784.             
  1785. if (kk<=k)
  1786.             { 
  1787. for (i=kk; i<=m; i++)
  1788.                 { 
  1789. ix=(i-1)*m+kk-1; 
  1790. iy=(i-1)*n+kk-1;
  1791.                     mtxU.m_pData[ix]=m_pData[iy];
  1792.                 }
  1793.             }
  1794.             
  1795. if (kk<=l)
  1796.             { 
  1797. d=0.0;
  1798.                 for (i=kk+1; i<=n; i++)
  1799. d=d+e[i-1]*e[i-1];
  1800.                 
  1801. e[kk-1]=sqrt(d);
  1802.                 if (e[kk-1]!=0.0)
  1803.                 { 
  1804. if (e[kk]!=0.0)
  1805.                     { 
  1806. e[kk-1]=fabs(e[kk-1]);
  1807.                         if (e[kk]<0.0) 
  1808. e[kk-1]=-e[kk-1];
  1809.                     }
  1811.                     for (i=kk+1; i<=n; i++)
  1812.                       e[i-1]=e[i-1]/e[kk-1];
  1813.                     
  1814. e[kk]=1.0+e[kk];
  1815.                 }
  1816.                 
  1817. e[kk-1]=-e[kk-1];
  1818.                 if ((kk+1<=m)&&(e[kk-1]!=0.0))
  1819.                 { 
  1820. for (i=kk+1; i<=m; i++) 
  1821. w[i-1]=0.0;
  1822.                     
  1823. for (j=kk+1; j<=n; j++)
  1824. for (i=kk+1; i<=m; i++)
  1825. w[i-1]=w[i-1]+e[j-1]*m_pData[(i-1)*n+j-1];
  1826.                     
  1827. for (j=kk+1; j<=n; j++)
  1828. {
  1829. for (i=kk+1; i<=m; i++)
  1830.                         { 
  1831. ix=(i-1)*n+j-1;
  1832. m_pData[ix]=m_pData[ix]-w[i-1]*e[j-1]/e[kk];
  1833.                         }
  1834. }
  1835.                 }
  1836.                 
  1837. for (i=kk+1; i<=n; i++)
  1838.                   mtxV.m_pData[(i-1)*n+kk-1]=e[i-1];
  1839.             }
  1840.         }
  1841.     }
  1842.     
  1843. mm=n;
  1844.     if (m+1<n) 
  1845. mm=m+1;
  1846.     if (k<n) 
  1847. s[k]=m_pData[k*n+k];
  1848.     if (m<mm) 
  1849. s[mm-1]=0.0;
  1850.     if (l+1<mm) 
  1851. e[l]=m_pData[l*n+mm-1];
  1853.     e[mm-1]=0.0;
  1854.     nn=m;
  1855.     if (m>n) 
  1856. nn=n;
  1857.     if (nn>=k+1)
  1858.     { 
  1859. for (j=k+1; j<=nn; j++)
  1860.         { 
  1861. for (i=1; i<=m; i++)
  1862. mtxU.m_pData[(i-1)*m+j-1]=0.0;
  1863.             mtxU.m_pData[(j-1)*m+j-1]=1.0;
  1864.         }
  1865.     }
  1866.     
  1867. if (k>=1)
  1868.     { 
  1869. for (ll=1; ll<=k; ll++)
  1870.         { 
  1871. kk=k-ll+1; 
  1872. iz=(kk-1)*m+kk-1;
  1873.             if (s[kk-1]!=0.0)
  1874.             { 
  1875. if (nn>=kk+1)
  1876. {
  1877. for (j=kk+1; j<=nn; j++)
  1878. { 
  1879. d=0.0;
  1880. for (i=kk; i<=m; i++)
  1881. { 
  1882. ix=(i-1)*m+kk-1;
  1883. iy=(i-1)*m+j-1;
  1884. d=d+mtxU.m_pData[ix]*mtxU.m_pData[iy]/mtxU.m_pData[iz];
  1885. }
  1887. d=-d;
  1888. for (i=kk; i<=m; i++)
  1889. { 
  1890. ix=(i-1)*m+j-1;
  1891. iy=(i-1)*m+kk-1;
  1892. mtxU.m_pData[ix]=mtxU.m_pData[ix]+d*mtxU.m_pData[iy];
  1893. }
  1894. }
  1895. }
  1896.                   
  1897. for (i=kk; i<=m; i++)
  1898. { 
  1899. ix=(i-1)*m+kk-1; 
  1900. mtxU.m_pData[ix]=-mtxU.m_pData[ix];
  1901. }
  1903. mtxU.m_pData[iz]=1.0+mtxU.m_pData[iz];
  1904. if (kk-1>=1)
  1905. {
  1906. for (i=1; i<=kk-1; i++)
  1907. mtxU.m_pData[(i-1)*m+kk-1]=0.0;
  1908. }
  1909. }
  1910.             else
  1911.             { 
  1912. for (i=1; i<=m; i++)
  1913. mtxU.m_pData[(i-1)*m+kk-1]=0.0;
  1914.                 mtxU.m_pData[(kk-1)*m+kk-1]=1.0;
  1915.             }
  1916. }
  1917.     }
  1919.     for (ll=1; ll<=n; ll++)
  1920.     { 
  1921. kk=n-ll+1; 
  1922. iz=kk*n+kk-1;
  1923.         
  1924. if ((kk<=l)&&(e[kk-1]!=0.0))
  1925.         { 
  1926. for (j=kk+1; j<=n; j++)
  1927.             { 
  1928. d=0.0;
  1929.                 for (i=kk+1; i<=n; i++)
  1930.                 { 
  1931. ix=(i-1)*n+kk-1; 
  1932. iy=(i-1)*n+j-1;
  1933.                     d=d+mtxV.m_pData[ix]*mtxV.m_pData[iy]/mtxV.m_pData[iz];
  1934.                 }
  1935.                 
  1936. d=-d;
  1937.                 for (i=kk+1; i<=n; i++)
  1938.                 { 
  1939. ix=(i-1)*n+j-1; 
  1940. iy=(i-1)*n+kk-1;
  1941.                     mtxV.m_pData[ix]=mtxV.m_pData[ix]+d*mtxV.m_pData[iy];
  1942.                 }
  1943.             }
  1944.         }
  1945.         
  1946. for (i=1; i<=n; i++)
  1947. mtxV.m_pData[(i-1)*n+kk-1]=0.0;
  1948.         
  1949. mtxV.m_pData[iz-n]=1.0;
  1950.     }
  1951.     
  1952. for (i=1; i<=m; i++)
  1953. for (j=1; j<=n; j++)
  1954. m_pData[(i-1)*n+j-1]=0.0;
  1955.     
  1956. m1=mm; 
  1957. it=60;
  1958.     while (TRUE)
  1959.     { 
  1960. if (mm==0)
  1961.         { 
  1962. ppp(m_pData,e,s,mtxV.m_pData,m,n);
  1963.             return TRUE;
  1964.         }
  1965.         if (it==0)
  1966.         { 
  1967. ppp(m_pData,e,s,mtxV.m_pData,m,n);
  1968.             return FALSE;
  1969.         }
  1970.         
  1971. kk=mm-1;
  1972. while ((kk!=0)&&(fabs(e[kk-1])!=0.0))
  1973.         { 
  1974. d=fabs(s[kk-1])+fabs(s[kk]);
  1975.             dd=fabs(e[kk-1]);
  1976.             if (dd>eps*d) 
  1977. kk=kk-1;
  1978.             else 
  1979. e[kk-1]=0.0;
  1980.         }
  1981.         
  1982. if (kk==mm-1)
  1983.         { 
  1984. kk=kk+1;
  1985.             if (s[kk-1]<0.0)
  1986.             { 
  1987. s[kk-1]=-s[kk-1];
  1988.                 for (i=1; i<=n; i++)
  1989.                 { 
  1990. ix=(i-1)*n+kk-1; 
  1991. mtxV.m_pData[ix]=-mtxV.m_pData[ix];}
  1992. }
  1994. while ((kk!=m1)&&(s[kk-1]<s[kk]))
  1995. { 
  1996. d=s[kk-1]; 
  1997. s[kk-1]=s[kk]; 
  1998. s[kk]=d;
  1999. if (kk<n)
  2000. {
  2001. for (i=1; i<=n; i++)
  2002. { 
  2003. ix=(i-1)*n+kk-1; 
  2004. iy=(i-1)*n+kk;
  2005. d=mtxV.m_pData[ix]; 
  2006. mtxV.m_pData[ix]=mtxV.m_pData[iy]; 
  2007. mtxV.m_pData[iy]=d;
  2008. }
  2009. }
  2011. if (kk<m)
  2012. {
  2013. for (i=1; i<=m; i++)
  2014. { 
  2015. ix=(i-1)*m+kk-1; 
  2016. iy=(i-1)*m+kk;
  2017. d=mtxU.m_pData[ix]; 
  2018. mtxU.m_pData[ix]=mtxU.m_pData[iy]; 
  2019. mtxU.m_pData[iy]=d;
  2020. }
  2021. }
  2023. kk=kk+1;
  2024.             }
  2025.             
  2026. it=60;
  2027.             mm=mm-1;
  2028.         }
  2029.         else
  2030.         { 
  2031. ks=mm;
  2032.             while ((ks>kk)&&(fabs(s[ks-1])!=0.0))
  2033.             { 
  2034. d=0.0;
  2035.                 if (ks!=mm) 
  2036. d=d+fabs(e[ks-1]);
  2037.                 if (ks!=kk+1) 
  2038. d=d+fabs(e[ks-2]);
  2039.                 
  2040. dd=fabs(s[ks-1]);
  2041.                 if (dd>eps*d) 
  2042. ks=ks-1;
  2043.                 else 
  2044. s[ks-1]=0.0;
  2045.             }
  2046.             
  2047. if (ks==kk)
  2048.             { 
  2049. kk=kk+1;
  2050.                 d=fabs(s[mm-1]);
  2051.                 t=fabs(s[mm-2]);
  2052.                 if (t>d) 
  2053. d=t;
  2054.                 
  2055. t=fabs(e[mm-2]);
  2056.                 if (t>d) 
  2057. d=t;
  2058.                 
  2059. t=fabs(s[kk-1]);
  2060.                 if (t>d) 
  2061. d=t;
  2062.                 
  2063. t=fabs(e[kk-1]);
  2064.                 if (t>d) 
  2065. d=t;
  2066.                 
  2067. sm=s[mm-1]/d; 
  2068. sm1=s[mm-2]/d;
  2069.                 em1=e[mm-2]/d;
  2070.                 sk=s[kk-1]/d; 
  2071. ek=e[kk-1]/d;
  2072.                 b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
  2073.                 c=sm*em1; 
  2074. c=c*c; 
  2075. shh=0.0;
  2077.                 if ((b!=0.0)||(c!=0.0))
  2078.                 { 
  2079. shh=sqrt(b*b+c);
  2080.                     if (b<0.0) 
  2081. shh=-shh;
  2083.                     shh=c/(b+shh);
  2084.                 }
  2085.                 
  2086. fg[0]=(sk+sm)*(sk-sm)-shh;
  2087.                 fg[1]=sk*ek;
  2088.                 for (i=kk; i<=mm-1; i++)
  2089.                 { 
  2090. sss(fg,cs);
  2091.                     if (i!=kk) 
  2092. e[i-2]=fg[0];
  2094.                     fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];
  2095.                     e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];
  2096.                     fg[1]=cs[1]*s[i];
  2097.                     s[i]=cs[0]*s[i];
  2099.                     if ((cs[0]!=1.0)||(cs[1]!=0.0))
  2100. {
  2101. for (j=1; j<=n; j++)
  2102.                         { 
  2103. ix=(j-1)*n+i-1;
  2104. iy=(j-1)*n+i;
  2105. d=cs[0]*mtxV.m_pData[ix]+cs[1]*mtxV.m_pData[iy];
  2106. mtxV.m_pData[iy]=-cs[1]*mtxV.m_pData[ix]+cs[0]*mtxV.m_pData[iy];
  2107. mtxV.m_pData[ix]=d;
  2108.                         }
  2109. }
  2111.                     sss(fg,cs);
  2112.                     s[i-1]=fg[0];
  2113.                     fg[0]=cs[0]*e[i-1]+cs[1]*s[i];
  2114.                     s[i]=-cs[1]*e[i-1]+cs[0]*s[i];
  2115.                     fg[1]=cs[1]*e[i];
  2116.                     e[i]=cs[0]*e[i];
  2118.                     if (i<m)
  2119. {
  2120. if ((cs[0]!=1.0)||(cs[1]!=0.0))
  2121. {
  2122. for (j=1; j<=m; j++)
  2123. { 
  2124. ix=(j-1)*m+i-1;
  2125. iy=(j-1)*m+i;
  2126. d=cs[0]*mtxU.m_pData[ix]+cs[1]*mtxU.m_pData[iy];
  2127. mtxU.m_pData[iy]=-cs[1]*mtxU.m_pData[ix]+cs[0]*mtxU.m_pData[iy];
  2128. mtxU.m_pData[ix]=d;
  2129. }
  2130. }
  2131. }
  2132.                 }
  2133.                 
  2134. e[mm-2]=fg[0];
  2135.                 it=it-1;
  2136.             }
  2137.             else
  2138.             { 
  2139. if (ks==mm)
  2140.                 { 
  2141. kk=kk+1;
  2142.                     fg[1]=e[mm-2]; 
  2143. e[mm-2]=0.0;
  2144.                     for (ll=kk; ll<=mm-1; ll++)
  2145.                     { 
  2146. i=mm+kk-ll-1;
  2147.                         fg[0]=s[i-1];
  2148.                         sss(fg,cs);
  2149.                         s[i-1]=fg[0];
  2150.                         if (i!=kk)
  2151.                         { 
  2152. fg[1]=-cs[1]*e[i-2];
  2153.                             e[i-2]=cs[0]*e[i-2];
  2154.                         }
  2155.                         
  2156. if ((cs[0]!=1.0)||(cs[1]!=0.0))
  2157. {
  2158. for (j=1; j<=n; j++)
  2159.                             { 
  2160. ix=(j-1)*n+i-1;
  2161. iy=(j-1)*n+mm-1;
  2162. d=cs[0]*mtxV.m_pData[ix]+cs[1]*mtxV.m_pData[iy];
  2163. mtxV.m_pData[iy]=-cs[1]*mtxV.m_pData[ix]+cs[0]*mtxV.m_pData[iy];
  2164. mtxV.m_pData[ix]=d;
  2165.                             }
  2166. }
  2167.                     }
  2168.                 }
  2169.                 else
  2170.                 { 
  2171. kk=ks+1;
  2172.                     fg[1]=e[kk-2];
  2173.                     e[kk-2]=0.0;
  2174.                     for (i=kk; i<=mm; i++)
  2175.                     { 
  2176. fg[0]=s[i-1];
  2177.                         sss(fg,cs);
  2178.                         s[i-1]=fg[0];
  2179.                         fg[1]=-cs[1]*e[i-1];
  2180.                         e[i-1]=cs[0]*e[i-1];
  2181.                         if ((cs[0]!=1.0)||(cs[1]!=0.0))
  2182. {
  2183. for (j=1; j<=m; j++)
  2184.                             { 
  2185. ix=(j-1)*m+i-1;
  2186. iy=(j-1)*m+kk-2;
  2187. d=cs[0]*mtxU.m_pData[ix]+cs[1]*mtxU.m_pData[iy];
  2188. mtxU.m_pData[iy]=-cs[1]*mtxU.m_pData[ix]+cs[0]*mtxU.m_pData[iy];
  2189. mtxU.m_pData[ix]=d;
  2190.                             }
  2191. }
  2192.                     }
  2193.                 }
  2194.             }
  2195.         }
  2196.     }
  2197.     
  2198. return TRUE;
  2199. }
  2201. //////////////////////////////////////////////////////////////////////
  2202. // 内部函数，由SplitUV函数调用
  2203. //////////////////////////////////////////////////////////////////////
  2204. void CMatrix::ppp(double a[], double e[], double s[], double v[], int m, int n)
  2205. { 
  2206. int i,j,p,q;
  2207.     double d;
  2209.     if (m>=n) 
  2210. i=n;
  2211.     else 
  2212. i=m;
  2214.     for (j=1; j<=i-1; j++)
  2215.     { 
  2216. a[(j-1)*n+j-1]=s[j-1];
  2217.         a[(j-1)*n+j]=e[j-1];
  2218.     }
  2219.     
  2220. a[(i-1)*n+i-1]=s[i-1];
  2221.     if (m<n) 
  2222. a[(i-1)*n+i]=e[i-1];
  2223.     
  2224. for (i=1; i<=n-1; i++)
  2225. {
  2226. for (j=i+1; j<=n; j++)
  2227. { 
  2228. p=(i-1)*n+j-1; 
  2229. q=(j-1)*n+i-1;
  2230. d=v[p]; 
  2231. v[p]=v[q]; 
  2232. v[q]=d;
  2233. }
  2234. }
  2235. }
  2237. //////////////////////////////////////////////////////////////////////
  2238. // 内部函数，由SplitUV函数调用
  2239. //////////////////////////////////////////////////////////////////////
  2240. void CMatrix::sss(double fg[2], double cs[2])
  2241. { 
  2242. double r,d;
  2243.     
  2244. if ((fabs(fg[0])+fabs(fg[1]))==0.0)
  2245.     { 
  2246. cs[0]=1.0; 
  2247. cs[1]=0.0; 
  2248. d=0.0;
  2249. }
  2250.     else 
  2251.     { 
  2252. d=sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
  2253.         if (fabs(fg[0])>fabs(fg[1]))
  2254.         { 
  2255. d=fabs(d);
  2256.             if (fg[0]<0.0) 
  2257. d=-d;
  2258.         }
  2259.         if (fabs(fg[1])>=fabs(fg[0]))
  2260.         { 
  2261. d=fabs(d);
  2262.             if (fg[1]<0.0) 
  2263. d=-d;
  2264.         }
  2265.         
  2266. cs[0]=fg[0]/d; 
  2267. cs[1]=fg[1]/d;
  2268.     }
  2269.     
  2270. r=1.0;
  2271.     if (fabs(fg[0])>fabs(fg[1])) 
  2272. r=cs[1];
  2273.     else if (cs[0]!=0.0) 
  2274. r=1.0/cs[0];
  2276.     fg[0]=d; 
  2277. fg[1]=r;
  2278. }
  2280. //////////////////////////////////////////////////////////////////////
  2281. // 求广义逆的奇异值分解法，分解成功后，原矩阵对角线元素就是矩阵的奇异值
  2282. //
  2283. // 参数：
  2284. // 1. CMatrix& mtxAP - 返回原矩阵的广义逆矩阵
  2285. // 2. CMatrix& mtxU - 返回分解后的U矩阵
  2286. // 3. CMatrix& mtxV - 返回分解后的V矩阵
  2287. // 4. double eps - 计算精度，默认值为0.000001
  2288. //
  2289. // 返回值：BOOL型，求解是否成功
  2290. //////////////////////////////////////////////////////////////////////
  2291. BOOL CMatrix::GInvertUV(CMatrix& mtxAP, CMatrix& mtxU, CMatrix& mtxV, double eps /*= 0.000001*/)
  2292. { 
  2293. int i,j,k,l,t,p,q,f;
  2295. // 调用奇异值分解
  2296.     if (! SplitUV(mtxU, mtxV, eps))
  2297. return FALSE;
  2299. int m = m_nNumRows;
  2300. int n = m_nNumColumns;
  2302. // 初始化广义逆矩阵
  2303. if (! mtxAP.Init(n, m))
  2304. return FALSE;
  2306. // 计算广义逆矩阵
  2308.     j=n;
  2309.     if (m<n) 
  2310. j=m;
  2311.     j=j-1;
  2312.     k=0;
  2313.     while ((k<=j)&&(m_pData[k*n+k]!=0.0)) 
  2314. k=k+1;
  2316.     k=k-1;
  2317.     for (i=0; i<=n-1; i++)
  2318. {
  2319. for (j=0; j<=m-1; j++)
  2320. { 
  2321. t=i*m+j;
  2322. mtxAP.m_pData[t]=0.0;
  2323. for (l=0; l<=k; l++)
  2324. { 
  2325. f=l*n+i; 
  2326. p=j*m+l; 
  2327. q=l*n+l;
  2328. mtxAP.m_pData[t]=mtxAP.m_pData[t]+mtxV.m_pData[f]*mtxU.m_pData[p]/m_pData[q];
  2329. }
  2330. }
  2331. }
  2333.     return TRUE;
  2334. }
  2336. //////////////////////////////////////////////////////////////////////
  2337. // 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法
  2338. //
  2339. // 参数：
  2340. // 1. CMatrix& mtxQ - 返回豪斯荷尔德变换的乘积矩阵Q
  2341. // 2. CMatrix& mtxT - 返回求得的对称三对角阵
  2342. // 3. double dblB[] - 一维数组，长度为矩阵的阶数，返回对称三对角阵的主对角线元素
  2343. // 4. double dblC[] - 一维数组，长度为矩阵的阶数，前n-1个元素返回对称三对角阵的次对角线元素
  2344. //
  2345. // 返回值：BOOL型，求解是否成功
  2346. //////////////////////////////////////////////////////////////////////
  2347. BOOL CMatrix::MakeSymTri(CMatrix& mtxQ, CMatrix& mtxT, double dblB[], double dblC[])
  2348. { 
  2349. int i,j,k,u;
  2350.     double h,f,g,h2;
  2351.     
  2352. // 初始化矩阵Q和T
  2353. if (! mtxQ.Init(m_nNumColumns, m_nNumColumns) ||
  2354. ! mtxT.Init(m_nNumColumns, m_nNumColumns))
  2355. return FALSE;
  2357. if (dblB == NULL || dblC == NULL)
  2358. return FALSE;
  2360. for (i=0; i<=m_nNumColumns-1; i++)
  2361. {
  2362. for (j=0; j<=m_nNumColumns-1; j++)
  2363. { 
  2364. u=i*m_nNumColumns+j; 
  2365. mtxQ.m_pData[u]=m_pData[u];
  2366. }
  2367. }
  2369.     for (i=m_nNumColumns-1; i>=1; i--)
  2370.     { 
  2371. h=0.0;
  2372.         if (i>1)
  2373. {
  2374. for (k=0; k<=i-1; k++)
  2375.             { 
  2376. u=i*m_nNumColumns+k; 
  2377. h=h+mtxQ.m_pData[u]*mtxQ.m_pData[u];
  2378. }
  2379. }
  2381.         if (h == 0.0)
  2382.         { 
  2383. dblC[i]=0.0;
  2384.             if (i==1) 
  2385. dblC[i]=mtxQ.m_pData[i*m_nNumColumns+i-1];
  2386.             dblB[i]=0.0;
  2387.         }
  2388.         else
  2389.         { 
  2390. dblC[i]=sqrt(h);
  2391.             u=i*m_nNumColumns+i-1;
  2392.             if (mtxQ.m_pData[u]>0.0) 
  2393. dblC[i]=-dblC[i];
  2395.             h=h-mtxQ.m_pData[u]*dblC[i];
  2396.             mtxQ.m_pData[u]=mtxQ.m_pData[u]-dblC[i];
  2397.             f=0.0;
  2398.             for (j=0; j<=i-1; j++)
  2399.             { 
  2400. mtxQ.m_pData[j*m_nNumColumns+i]=mtxQ.m_pData[i*m_nNumColumns+j]/h;
  2401.                 g=0.0;
  2402.                 for (k=0; k<=j; k++)
  2403. g=g+mtxQ.m_pData[j*m_nNumColumns+k]*mtxQ.m_pData[i*m_nNumColumns+k];
  2405. if (j+1<=i-1)
  2406. for (k=j+1; k<=i-1; k++)
  2407. g=g+mtxQ.m_pData[k*m_nNumColumns+j]*mtxQ.m_pData[i*m_nNumColumns+k];
  2409.                 dblC[j]=g/h;
  2410.                 f=f+g*mtxQ.m_pData[j*m_nNumColumns+i];
  2411.             }
  2412.             
  2413. h2=f/(h+h);
  2414.             for (j=0; j<=i-1; j++)
  2415.             { 
  2416. f=mtxQ.m_pData[i*m_nNumColumns+j];
  2417.                 g=dblC[j]-h2*f;
  2418.                 dblC[j]=g;
  2419.                 for (k=0; k<=j; k++)
  2420.                 { 
  2421. u=j*m_nNumColumns+k;
  2422.                     mtxQ.m_pData[u]=mtxQ.m_pData[u]-f*dblC[k]-g*mtxQ.m_pData[i*m_nNumColumns+k];
  2423.                 }
  2424.             }
  2425.             
  2426. dblB[i]=h;
  2427.         }
  2428.     }
  2429.     
  2430. for (i=0; i<=m_nNumColumns-2; i++) 
  2431. dblC[i]=dblC[i+1];
  2432.     
  2433. dblC[m_nNumColumns-1]=0.0;
  2434.     dblB[0]=0.0;
  2435.     for (i=0; i<=m_nNumColumns-1; i++)
  2436.     { 
  2437. if ((dblB[i]!=0.0)&&(i-1>=0))
  2438. {
  2439. for (j=0; j<=i-1; j++)
  2440.             { 
  2441. g=0.0;
  2442. for (k=0; k<=i-1; k++)
  2443. g=g+mtxQ.m_pData[i*m_nNumColumns+k]*mtxQ.m_pData[k*m_nNumColumns+j];
  2445. for (k=0; k<=i-1; k++)
  2446.                 { 
  2447. u=k*m_nNumColumns+j;
  2448. mtxQ.m_pData[u]=mtxQ.m_pData[u]-g*mtxQ.m_pData[k*m_nNumColumns+i];
  2449.                 }
  2450.             }
  2451. }
  2453.         u=i*m_nNumColumns+i;
  2454.         dblB[i]=mtxQ.m_pData[u]; mtxQ.m_pData[u]=1.0;
  2455.         if (i-1>=0)
  2456. {
  2457. for (j=0; j<=i-1; j++)
  2458.             { 
  2459. mtxQ.m_pData[i*m_nNumColumns+j]=0.0; 
  2460. mtxQ.m_pData[j*m_nNumColumns+i]=0.0;
  2461. }
  2462. }
  2463.     }
  2465.     // 构造对称三对角矩阵
  2466.     for (i=0; i<m_nNumColumns; ++i)
  2467. {
  2468.     for (j=0; j<m_nNumColumns; ++j)
  2469. {
  2470.             mtxT.SetElement(i, j, 0);
  2471.             k = i - j;
  2472.             if (k == 0) 
  2473.             mtxT.SetElement(i, j, dblB[j]);
  2474. else if (k == 1)
  2475.             mtxT.SetElement(i, j, dblC[j]);
  2476. else if (k == -1)
  2477.             mtxT.SetElement(i, j, dblC[i]);
  2478.         }
  2479.     }
  2481. return TRUE;
  2482. }
  2484. //////////////////////////////////////////////////////////////////////
  2485. // 实对称三对角阵的全部特征值与特征向量的计算
  2486. //
  2487. // 参数：
  2488. // 1. double dblB[] - 一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；
  2489. //    返回时存放全部特征值。
  2490. // 2. double dblC[] - 一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的次对角线元素
  2491. // 3. CMatrix& mtxQ - 如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；
  2492. //    如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积矩阵Q，则返回矩阵A的
  2493. //    特征值向量矩阵。其中第i列为与数组dblB中第j个特征值对应的特征向量。
  2494. // 4. int nMaxIt - 迭代次数，默认值为60
  2495. // 5. double eps - 计算精度，默认值为0.000001
  2496. //
  2497. // 返回值：BOOL型，求解是否成功
  2498. //////////////////////////////////////////////////////////////////////
  2499. BOOL CMatrix::SymTriEigenv(double dblB[], double dblC[], CMatrix& mtxQ, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  2500. {
  2501. int i,j,k,m,it,u,v;
  2502.     double d,f,h,g,p,r,e,s;
  2503.     
  2504. // 初值
  2505. int n = mtxQ.GetNumColumns();
  2506. dblC[n-1]=0.0; 
  2507. d=0.0; 
  2508. f=0.0;
  2509.     
  2510. // 迭代计算
  2512. for (j=0; j<=n-1; j++)
  2513.     { 
  2514. it=0;
  2515.         h=eps*(fabs(dblB[j])+fabs(dblC[j]));
  2516.         if (h>d) 
  2517. d=h;
  2518.         
  2519. m=j;
  2520.         while ((m<=n-1)&&(fabs(dblC[m])>d)) 
  2521. m=m+1;
  2522.         
  2523. if (m!=j)
  2524.         { 
  2525. do
  2526.             { 
  2527. if (it==nMaxIt)
  2528. return FALSE;
  2530.                 it=it+1;
  2531.                 g=dblB[j];
  2532.                 p=(dblB[j+1]-g)/(2.0*dblC[j]);
  2533.                 r=sqrt(p*p+1.0);
  2534.                 if (p>=0.0) 
  2535. dblB[j]=dblC[j]/(p+r);
  2536.                 else 
  2537. dblB[j]=dblC[j]/(p-r);
  2538.                 
  2539. h=g-dblB[j];
  2540.                 for (i=j+1; i<=n-1; i++)
  2541. dblB[i]=dblB[i]-h;
  2542.                 
  2543. f=f+h; 
  2544. p=dblB[m]; 
  2545. e=1.0; 
  2546. s=0.0;
  2547.                 for (i=m-1; i>=j; i--)
  2548.                 { 
  2549. g=e*dblC[i]; 
  2550. h=e*p;
  2551.                     if (fabs(p)>=fabs(dblC[i]))
  2552.                     { 
  2553. e=dblC[i]/p; 
  2554. r=sqrt(e*e+1.0);
  2555.                         dblC[i+1]=s*p*r; 
  2556. s=e/r; 
  2557. e=1.0/r;
  2558.                     }
  2559.                     else
  2560. { 
  2561. e=p/dblC[i]; 
  2562. r=sqrt(e*e+1.0);
  2563.                         dblC[i+1]=s*dblC[i]*r;
  2564.                         s=1.0/r; 
  2565. e=e/r;
  2566.                     }
  2567.                     
  2568. p=e*dblB[i]-s*g;
  2569.                     dblB[i+1]=h+s*(e*g+s*dblB[i]);
  2570.                     for (k=0; k<=n-1; k++)
  2571.                     { 
  2572. u=k*n+i+1; 
  2573. v=u-1;
  2574.                         h=mtxQ.m_pData[u]; 
  2575. mtxQ.m_pData[u]=s*mtxQ.m_pData[v]+e*h;
  2576.                         mtxQ.m_pData[v]=e*mtxQ.m_pData[v]-s*h;
  2577.                     }
  2578.                 }
  2579.                 
  2580. dblC[j]=s*p; 
  2581. dblB[j]=e*p;
  2582.             
  2583. } while (fabs(dblC[j])>d);
  2584.         }
  2585.         
  2586. dblB[j]=dblB[j]+f;
  2587.     }
  2588.     
  2589. for (i=0; i<=n-1; i++)
  2590.     { 
  2591. k=i; 
  2592. p=dblB[i];
  2593.         if (i+1<=n-1)
  2594.         { 
  2595. j=i+1;
  2596.             while ((j<=n-1)&&(dblB[j]<=p))
  2597.             { 
  2598. k=j; 
  2599. p=dblB[j]; 
  2600. j=j+1;
  2601. }
  2602.         }
  2604.         if (k!=i)
  2605.         { 
  2606. dblB[k]=dblB[i]; 
  2607. dblB[i]=p;
  2608.             for (j=0; j<=n-1; j++)
  2609.             { 
  2610. u=j*n+i; 
  2611. v=j*n+k;
  2612.                 p=mtxQ.m_pData[u]; 
  2613. mtxQ.m_pData[u]=mtxQ.m_pData[v]; 
  2614. mtxQ.m_pData[v]=p;
  2615.             }
  2616.         }
  2617.     }
  2618.     
  2619. return TRUE;
  2620. }
  2622. //////////////////////////////////////////////////////////////////////
  2623. // 约化一般实矩阵为赫申伯格矩阵的初等相似变换法
  2624. //
  2625. // 参数：无
  2626. //
  2627. // 返回值：无
  2628. //////////////////////////////////////////////////////////////////////
  2629. void CMatrix::MakeHberg()
  2630. { 
  2631. int i,j,k,u,v;
  2632.     double d,t;
  2634.     for (k=1; k<=m_nNumColumns-2; k++)
  2635.     { 
  2636. d=0.0;
  2637.         for (j=k; j<=m_nNumColumns-1; j++)
  2638.         { 
  2639. u=j*m_nNumColumns+k-1; 
  2640. t=m_pData[u];
  2641.             if (fabs(t)>fabs(d))
  2642.             { 
  2643. d=t; 
  2644. i=j;
  2645. }
  2646.         }
  2647.         
  2648. if (d != 0.0)
  2649.         { 
  2650. if (i!=k)
  2651.             { 
  2652. for (j=k-1; j<=m_nNumColumns-1; j++)
  2653.                 { 
  2654. u=i*m_nNumColumns+j; 
  2655. v=k*m_nNumColumns+j;
  2656.                     t=m_pData[u]; 
  2657. m_pData[u]=m_pData[v]; 
  2658. m_pData[v]=t;
  2659.                 }
  2660.                 
  2661. for (j=0; j<=m_nNumColumns-1; j++)
  2662.                 { 
  2663. u=j*m_nNumColumns+i; 
  2664. v=j*m_nNumColumns+k;
  2665.                     t=m_pData[u]; 
  2666. m_pData[u]=m_pData[v]; 
  2667. m_pData[v]=t;
  2668.                 }
  2669.             }
  2670.             
  2671. for (i=k+1; i<=m_nNumColumns-1; i++)
  2672.             { 
  2673. u=i*m_nNumColumns+k-1; 
  2674. t=m_pData[u]/d; 
  2675. m_pData[u]=0.0;
  2676.                 for (j=k; j<=m_nNumColumns-1; j++)
  2677.                 { 
  2678. v=i*m_nNumColumns+j;
  2679.                     m_pData[v]=m_pData[v]-t*m_pData[k*m_nNumColumns+j];
  2680.                 }
  2681.                 
  2682. for (j=0; j<=m_nNumColumns-1; j++)
  2683.                 { 
  2684. v=j*m_nNumColumns+k;
  2685.                     m_pData[v]=m_pData[v]+t*m_pData[j*m_nNumColumns+i];
  2686.                 }
  2687.             }
  2688.         }
  2689.     }
  2690. }
  2692. //////////////////////////////////////////////////////////////////////
  2693. // 求赫申伯格矩阵全部特征值的QR方法
  2694. //
  2695. // 参数：
  2696. // 1. double dblU[] - 一维数组，长度为矩阵的阶数，返回时存放特征值的实部
  2697. // 2. double dblV[] - 一维数组，长度为矩阵的阶数，返回时存放特征值的虚部
  2698. // 3. int nMaxIt - 迭代次数，默认值为60
  2699. // 4. double eps - 计算精度，默认值为0.000001
  2700. //
  2701. // 返回值：BOOL型，求解是否成功
  2702. //////////////////////////////////////////////////////////////////////
  2703. BOOL CMatrix::HBergEigenv(double dblU[], double dblV[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  2704. { 
  2705. int m,it,i,j,k,l,ii,jj,kk,ll;
  2706.     double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
  2707.     
  2708. int n = m_nNumColumns;
  2710. it=0; 
  2711. m=n;
  2712.     while (m!=0)
  2713.     { 
  2714. l=m-1;
  2715.         while ((l>0)&&(fabs(m_pData[l*n+l-1]) > 
  2716. eps*(fabs(m_pData[(l-1)*n+l-1])+fabs(m_pData[l*n+l])))) 
  2717.   l=l-1;
  2719.         ii=(m-1)*n+m-1; 
  2720. jj=(m-1)*n+m-2;
  2721.         kk=(m-2)*n+m-1; 
  2722. ll=(m-2)*n+m-2;
  2723.         if (l==m-1)
  2724.         { 
  2725. dblU[m-1]=m_pData[(m-1)*n+m-1]; 
  2726. dblV[m-1]=0.0;
  2727.             m=m-1; 
  2728. it=0;
  2729.         }
  2730.         else if (l==m-2)
  2731.         { 
  2732. b=-(m_pData[ii]+m_pData[ll]);
  2733.             c=m_pData[ii]*m_pData[ll]-m_pData[jj]*m_pData[kk];
  2734.             w=b*b-4.0*c;
  2735.             y=sqrt(fabs(w));
  2736.             if (w>0.0)
  2737.             { 
  2738. xy=1.0;
  2739.                 if (b<0.0) 
  2740. xy=-1.0;
  2741.                 dblU[m-1]=(-b-xy*y)/2.0;
  2742.                 dblU[m-2]=c/dblU[m-1];
  2743.                 dblV[m-1]=0.0; dblV[m-2]=0.0;
  2744.             }
  2745.             else
  2746.             { 
  2747. dblU[m-1]=-b/2.0; 
  2748. dblU[m-2]=dblU[m-1];
  2749.                 dblV[m-1]=y/2.0; 
  2750. dblV[m-2]=-dblV[m-1];
  2751.             }
  2752.             
  2753. m=m-2; 
  2754. it=0;
  2755.         }
  2756.         else
  2757.         { 
  2758. if (it>=nMaxIt)
  2759. return FALSE;
  2761.             it=it+1;
  2762.             for (j=l+2; j<=m-1; j++)
  2763. m_pData[j*n+j-2]=0.0;
  2764.             for (j=l+3; j<=m-1; j++)
  2765. m_pData[j*n+j-3]=0.0;
  2766.             for (k=l; k<=m-2; k++)
  2767.             { 
  2768. if (k!=l)
  2769.                 { 
  2770. p=m_pData[k*n+k-1]; 
  2771. q=m_pData[(k+1)*n+k-1];
  2772.                     r=0.0;
  2773.                     if (k!=m-2) 
  2774. r=m_pData[(k+2)*n+k-1];
  2775.                 }
  2776.                 else
  2777.                 { 
  2778. x=m_pData[ii]+m_pData[ll];
  2779.                     y=m_pData[ll]*m_pData[ii]-m_pData[kk]*m_pData[jj];
  2780.                     ii=l*n+l; 
  2781. jj=l*n+l+1;
  2782.                     kk=(l+1)*n+l; 
  2783. ll=(l+1)*n+l+1;
  2784.                     p=m_pData[ii]*(m_pData[ii]-x)+m_pData[jj]*m_pData[kk]+y;
  2785.                     q=m_pData[kk]*(m_pData[ii]+m_pData[ll]-x);
  2786.                     r=m_pData[kk]*m_pData[(l+2)*n+l+1];
  2787.                 }
  2788.                 
  2789. if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
  2790.                 { 
  2791. xy=1.0;
  2792.                     if (p<0.0) 
  2793. xy=-1.0;
  2794.                     s=xy*sqrt(p*p+q*q+r*r);
  2795.                     if (k!=l) 
  2796. m_pData[k*n+k-1]=-s;
  2797.                     e=-q/s; 
  2798. f=-r/s; 
  2799. x=-p/s;
  2800.                     y=-x-f*r/(p+s);
  2801.                     g=e*r/(p+s);
  2802.                     z=-x-e*q/(p+s);
  2803.                     for (j=k; j<=m-1; j++)
  2804.                     { 
  2805. ii=k*n+j; 
  2806. jj=(k+1)*n+j;
  2807.                         p=x*m_pData[ii]+e*m_pData[jj];
  2808.                         q=e*m_pData[ii]+y*m_pData[jj];
  2809.                         r=f*m_pData[ii]+g*m_pData[jj];
  2810.                         if (k!=m-2)
  2811.                         { 
  2812. kk=(k+2)*n+j;
  2813.                             p=p+f*m_pData[kk];
  2814.                             q=q+g*m_pData[kk];
  2815.                             r=r+z*m_pData[kk]; 
  2816. m_pData[kk]=r;
  2817.                         }
  2818.                         
  2819. m_pData[jj]=q; m_pData[ii]=p;
  2820.                     }
  2821.                     
  2822. j=k+3;
  2823.                     if (j>=m-1) 
  2824. j=m-1;
  2825.                     
  2826. for (i=l; i<=j; i++)
  2827.                     { 
  2828. ii=i*n+k; 
  2829. jj=i*n+k+1;
  2830.                         p=x*m_pData[ii]+e*m_pData[jj];
  2831.                         q=e*m_pData[ii]+y*m_pData[jj];
  2832.                         r=f*m_pData[ii]+g*m_pData[jj];
  2833.                         if (k!=m-2)
  2834.                         { 
  2835. kk=i*n+k+2;
  2836.                             p=p+f*m_pData[kk];
  2837.                             q=q+g*m_pData[kk];
  2838.                             r=r+z*m_pData[kk]; 
  2839. m_pData[kk]=r;
  2840.                         }
  2841.                         
  2842. m_pData[jj]=q; 
  2843. m_pData[ii]=p;
  2844.                     }
  2845.                 }
  2846.             }
  2847.         }
  2848.     }
  2849.     
  2850. return TRUE;
  2851. }
  2853. //////////////////////////////////////////////////////////////////////
  2854. // 求实对称矩阵特征值与特征向量的雅可比法
  2855. //
  2856. // 参数：
  2857. // 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值
  2858. // 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与
  2859. //    数组dblEigenValue中第j个特征值对应的特征向量
  2860. // 3. int nMaxIt - 迭代次数，默认值为60
  2861. // 4. double eps - 计算精度，默认值为0.000001
  2862. //
  2863. // 返回值：BOOL型，求解是否成功
  2864. //////////////////////////////////////////////////////////////////////
  2865. BOOL CMatrix::JacobiEigenv(double dblEigenValue[], CMatrix& mtxEigenVector, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  2866. { 
  2867. int i,j,p,q,u,w,t,s,l;
  2868.     double fm,cn,sn,omega,x,y,d;
  2869.     
  2870. if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
  2871. return FALSE;
  2873. l=1;
  2874.     for (i=0; i<=m_nNumColumns-1; i++)
  2875.     { 
  2876. mtxEigenVector.m_pData[i*m_nNumColumns+i]=1.0;
  2877.         for (j=0; j<=m_nNumColumns-1; j++)
  2878. if (i!=j) 
  2879. mtxEigenVector.m_pData[i*m_nNumColumns+j]=0.0;
  2880.     }
  2881.     
  2882. while (TRUE)
  2883.     { 
  2884. fm=0.0;
  2885.         for (i=1; i<=m_nNumColumns-1; i++)
  2886. {
  2887. for (j=0; j<=i-1; j++)
  2888. { 
  2889. d=fabs(m_pData[i*m_nNumColumns+j]);
  2890. if ((i!=j)&&(d>fm))
  2891. { 
  2892. fm=d; 
  2893. p=i; 
  2894. q=j;
  2895. }
  2896. }
  2897. }
  2899.         if (fm<eps)
  2900. {
  2901. for (i=0; i<m_nNumColumns; ++i)
  2902. dblEigenValue[i] = GetElement(i,i);
  2903. return TRUE;
  2904. }
  2906.         if (l>nMaxIt)  
  2907. return FALSE;
  2908.         
  2909. l=l+1;
  2910.         u=p*m_nNumColumns+q; 
  2911. w=p*m_nNumColumns+p; 
  2912. t=q*m_nNumColumns+p; 
  2913. s=q*m_nNumColumns+q;
  2914.         x=-m_pData[u]; 
  2915. y=(m_pData[s]-m_pData[w])/2.0;
  2916.         omega=x/sqrt(x*x+y*y);
  2918.         if (y<0.0) 
  2919. omega=-omega;
  2921.         sn=1.0+sqrt(1.0-omega*omega);
  2922.         sn=omega/sqrt(2.0*sn);
  2923.         cn=sqrt(1.0-sn*sn);
  2924.         fm=m_pData[w];
  2925.         m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData[u]*omega;
  2926.         m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData[u]*omega;
  2927.         m_pData[u]=0.0; 
  2928. m_pData[t]=0.0;
  2929.         for (j=0; j<=m_nNumColumns-1; j++)
  2930. {
  2931. if ((j!=p)&&(j!=q))
  2932. { 
  2933. u=p*m_nNumColumns+j; w=q*m_nNumColumns+j;
  2934. fm=m_pData[u];
  2935. m_pData[u]=fm*cn+m_pData[w]*sn;
  2936. m_pData[w]=-fm*sn+m_pData[w]*cn;
  2937. }
  2938. }
  2940.         for (i=0; i<=m_nNumColumns-1; i++)
  2941. {
  2942. if ((i!=p)&&(i!=q))
  2943.             { 
  2944. u=i*m_nNumColumns+p; 
  2945. w=i*m_nNumColumns+q;
  2946. fm=m_pData[u];
  2947. m_pData[u]=fm*cn+m_pData[w]*sn;
  2948. m_pData[w]=-fm*sn+m_pData[w]*cn;
  2949.             }
  2950. }
  2952.         for (i=0; i<=m_nNumColumns-1; i++)
  2953.         { 
  2954. u=i*m_nNumColumns+p; 
  2955. w=i*m_nNumColumns+q;
  2956.             fm=mtxEigenVector.m_pData[u];
  2957.             mtxEigenVector.m_pData[u]=fm*cn+mtxEigenVector.m_pData[w]*sn;
  2958.             mtxEigenVector.m_pData[w]=-fm*sn+mtxEigenVector.m_pData[w]*cn;
  2959.         }
  2960.     }
  2961.     
  2962. for (i=0; i<m_nNumColumns; ++i)
  2963. dblEigenValue[i] = GetElement(i,i);
  2965. return TRUE;
  2966. }
  2968. //////////////////////////////////////////////////////////////////////
  2969. // 求实对称矩阵特征值与特征向量的雅可比过关法
  2970. //
  2971. // 参数：
  2972. // 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值
  2973. // 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与
  2974. //    数组dblEigenValue中第j个特征值对应的特征向量
  2975. // 3. double eps - 计算精度，默认值为0.000001
  2976. //
  2977. // 返回值：BOOL型，求解是否成功
  2978. //////////////////////////////////////////////////////////////////////
  2979. BOOL CMatrix::JacobiEigenv2(double dblEigenValue[], CMatrix& mtxEigenVector, double eps /*= 0.000001*/)
  2980. { 
  2981. int i,j,p,q,u,w,t,s;
  2982.     double ff,fm,cn,sn,omega,x,y,d;
  2983.     
  2984. if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
  2985. return FALSE;
  2987. for (i=0; i<=m_nNumColumns-1; i++)
  2988.     { 
  2989. mtxEigenVector.m_pData[i*m_nNumColumns+i]=1.0;
  2990.         for (j=0; j<=m_nNumColumns-1; j++)
  2991. if (i!=j) 
  2992. mtxEigenVector.m_pData[i*m_nNumColumns+j]=0.0;
  2993.     }
  2994.     
  2995. ff=0.0;
  2996.     for (i=1; i<=m_nNumColumns-1; i++)
  2997. {
  2998. for (j=0; j<=i-1; j++)
  2999. { 
  3000. d=m_pData[i*m_nNumColumns+j]; 
  3001. ff=ff+d*d; 
  3002. }
  3003. }
  3005.     ff=sqrt(2.0*ff);
  3007. Loop_0:
  3008.     
  3009. ff=ff/(1.0*m_nNumColumns);
  3011. Loop_1:
  3013.     for (i=1; i<=m_nNumColumns-1; i++)
  3014. {
  3015. for (j=0; j<=i-1; j++)
  3016.         { 
  3017. d=fabs(m_pData[i*m_nNumColumns+j]);
  3018.             if (d>ff)
  3019.             { 
  3020. p=i; 
  3021. q=j;
  3022.                 goto Loop_2;
  3023.             }
  3024.         }
  3025. }
  3026.         
  3027. if (ff<eps) 
  3028. {
  3029. for (i=0; i<m_nNumColumns; ++i)
  3030. dblEigenValue[i] = GetElement(i,i);
  3031. return TRUE;
  3032. }
  3033.     
  3034. goto Loop_0;
  3036. Loop_2: 
  3038. u=p*m_nNumColumns+q; 
  3039. w=p*m_nNumColumns+p; 
  3040. t=q*m_nNumColumns+p; 
  3041. s=q*m_nNumColumns+q;
  3042.     x=-m_pData[u]; 
  3043. y=(m_pData[s]-m_pData[w])/2.0;
  3044.     omega=x/sqrt(x*x+y*y);
  3045.     if (y<0.0) 
  3046. omega=-omega;
  3047.     
  3048. sn=1.0+sqrt(1.0-omega*omega);
  3049.     sn=omega/sqrt(2.0*sn);
  3050.     cn=sqrt(1.0-sn*sn);
  3051.     fm=m_pData[w];
  3052.     m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData[u]*omega;
  3053.     m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData[u]*omega;
  3054.     m_pData[u]=0.0; m_pData[t]=0.0;
  3055.     
  3056. for (j=0; j<=m_nNumColumns-1; j++)
  3057. {
  3058. if ((j!=p)&&(j!=q))
  3059. { 
  3060. u=p*m_nNumColumns+j; 
  3061. w=q*m_nNumColumns+j;
  3062. fm=m_pData[u];
  3063. m_pData[u]=fm*cn+m_pData[w]*sn;
  3064. m_pData[w]=-fm*sn+m_pData[w]*cn;
  3065. }
  3066. }
  3068.     for (i=0; i<=m_nNumColumns-1; i++)
  3069.     {
  3070. if ((i!=p)&&(i!=q))
  3071.         { 
  3072. u=i*m_nNumColumns+p; 
  3073. w=i*m_nNumColumns+q;
  3074. fm=m_pData[u];
  3075. m_pData[u]=fm*cn+m_pData[w]*sn;
  3076. m_pData[w]=-fm*sn+m_pData[w]*cn;
  3077.         }
  3078. }
  3079.     
  3080. for (i=0; i<=m_nNumColumns-1; i++)
  3081.     { 
  3082. u=i*m_nNumColumns+p; 
  3083. w=i*m_nNumColumns+q;
  3084.         fm=mtxEigenVector.m_pData[u];
  3085.         mtxEigenVector.m_pData[u]=fm*cn+mtxEigenVector.m_pData[w]*sn;
  3086.         mtxEigenVector.m_pData[w]=-fm*sn+mtxEigenVector.m_pData[w]*cn;
  3087. }
  3089. goto Loop_1;
  3090. }