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

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

数学计算

开发平台：

Visual C++

LEquations.cpp：源码内容
  1. //////////////////////////////////////////////////////////////////////
  2. // LEquations.cpp
  3. //
  4. // 求解线性方程组的类 CLEquations 的实现代码
  5. //
  6. // 周长发编制, 2002/8
  7. //////////////////////////////////////////////////////////////////////
  9. #include "stdafx.h"
  10. #include "LEquations.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. CLEquations::CLEquations()
  26. {
  27. }
  29. //////////////////////////////////////////////////////////////////////
  30. // 指定系数和常数构造函数
  31. //
  32. // 参数：
  33. // 1. const CMatrix& mtxCoef - 指定的系数矩阵
  34. // 2. const CMatrix& mtxConst - 指定的常数矩阵
  35. //////////////////////////////////////////////////////////////////////
  36. CLEquations::CLEquations(const CMatrix& mtxCoef, const CMatrix& mtxConst)
  37. {
  38. ASSERT(Init(mtxCoef, mtxConst));
  39. }
  41. //////////////////////////////////////////////////////////////////////
  42. // 析构函数
  43. //////////////////////////////////////////////////////////////////////
  44. CLEquations::~CLEquations()
  45. {
  46. }
  48. //////////////////////////////////////////////////////////////////////
  49. // 初始化函数
  50. //
  51. // 参数：
  52. // 1. const CMatrix& mtxCoef - 指定的系数矩阵
  53. // 2. const CMatrix& mtxConst - 指定的常数矩阵
  54. //
  55. // 返回值：BOOL 型，初始化是否成功
  56. //////////////////////////////////////////////////////////////////////
  57. BOOL CLEquations::Init(const CMatrix& mtxCoef, const CMatrix& mtxConst)
  58. {
  59. if (mtxCoef.GetNumRows() != mtxConst.GetNumRows())
  60. return FALSE;
  62. m_mtxCoef = mtxCoef;
  63. m_mtxConst = mtxConst;
  65. return TRUE;
  66. }
  68. //////////////////////////////////////////////////////////////////////
  69. // 获取系数矩阵
  70. //
  71. // 参数：无
  72. //
  73. // 返回值：CMatrix 型，返回系数矩阵
  74. //////////////////////////////////////////////////////////////////////
  75. inline CMatrix CLEquations::GetCoefMatrix() const
  76. {
  77. return m_mtxCoef;
  78. }
  80. //////////////////////////////////////////////////////////////////////
  81. // 获取常数矩阵
  82. //
  83. // 参数：无
  84. //
  85. // 返回值：CMatrix 型，返回常数矩阵
  86. //////////////////////////////////////////////////////////////////////
  87. inline CMatrix CLEquations::GetConstMatrix() const
  88. {
  89. return m_mtxConst;
  90. }
  92. //////////////////////////////////////////////////////////////////////
  93. // 获取方程个数
  94. //
  95. // 参数：无
  96. //
  97. // 返回值：int 型，返回方程组方程的个数
  98. //////////////////////////////////////////////////////////////////////
  99. inline int CLEquations::GetNumEquations() const
  100. {
  101. return GetCoefMatrix().GetNumRows();
  102. }
  104. //////////////////////////////////////////////////////////////////////
  105. // 获取未知数个数
  106. //
  107. // 参数：无
  108. //
  109. // 返回值：int 型，返回方程组未知数的个数
  110. //////////////////////////////////////////////////////////////////////
  111. inline int CLEquations::GetNumUnknowns() const
  112. {
  113. return GetCoefMatrix().GetNumColumns();
  114. }
  116. //////////////////////////////////////////////////////////////////////
  117. // 全选主元高斯消去法
  118. //
  119. // 参数：
  120. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组的解
  121. //
  122. // 返回值：BOOL 型，方程组求解是否成功
  123. //////////////////////////////////////////////////////////////////////
  124. BOOL CLEquations::GetRootsetGauss(CMatrix& mtxResult)
  125. { 
  126. int *pnJs,l,k,i,j,nIs,p,q;
  127.     double d,t;
  129. // 方程组的属性，将常数矩阵赋给解矩阵
  130. mtxResult = m_mtxConst;
  131. double *pDataCoef = m_mtxCoef.GetData();
  132. double *pDataConst = mtxResult.GetData();
  133. int n = GetNumUnknowns();
  135. // 临时缓冲区，存放列数
  136.     pnJs = new int[n];
  138. // 消元
  139.     l=1;
  140.     for (k=0;k<=n-2;k++)
  141.     { 
  142. d=0.0;
  143.         for (i=k;i<=n-1;i++)
  144. {
  145. for (j=k;j<=n-1;j++)
  146.             { 
  147. t=fabs(pDataCoef[i*n+j]);
  148. if (t>d) 
  149. { 
  150. d=t; 
  151. pnJs[k]=j; 
  152. nIs=i;
  153. }
  154.             }
  155. }
  157.         if (d == 0.0) 
  158. l=0;
  159.         else
  160.         { 
  161. if (pnJs[k]!=k)
  162.             {
  163. for (i=0;i<=n-1;i++)
  164.                 { 
  165. p=i*n+k; 
  166. q=i*n+pnJs[k];
  167. t=pDataCoef[p]; 
  168. pDataCoef[p]=pDataCoef[q]; 
  169. pDataCoef[q]=t;
  170.                 }
  171. }
  173.             if (nIs!=k)
  174.             { 
  175. for (j=k;j<=n-1;j++)
  176.                 { 
  177. p=k*n+j; 
  178. q=nIs*n+j;
  179.                     t=pDataCoef[p]; 
  180. pDataCoef[p]=pDataCoef[q]; 
  181. pDataCoef[q]=t;
  182.                 }
  183.                 
  184. t=pDataConst[k]; 
  185. pDataConst[k]=pDataConst[nIs]; 
  186. pDataConst[nIs]=t;
  187.             }
  188.         }
  189.         
  190. // 求解失败
  191. if (l==0)
  192.         { 
  193. delete[] pnJs;
  194.             return FALSE;
  195.         }
  196.         
  197. d=pDataCoef[k*n+k];
  198.         for (j=k+1;j<=n-1;j++)
  199.         { 
  200. p=k*n+j; 
  201. pDataCoef[p]=pDataCoef[p]/d;
  202. }
  203.         
  204. pDataConst[k]=pDataConst[k]/d;
  205.         for (i=k+1;i<=n-1;i++)
  206.         { 
  207. for (j=k+1;j<=n-1;j++)
  208.             { 
  209. p=i*n+j;
  210.                 pDataCoef[p]=pDataCoef[p]-pDataCoef[i*n+k]*pDataCoef[k*n+j];
  211.             }
  212.             
  213. pDataConst[i]=pDataConst[i]-pDataCoef[i*n+k]*pDataConst[k];
  214.         }
  215.     }
  216.     
  217. // 求解失败
  218. d=pDataCoef[(n-1)*n+n-1];
  219.     if (d == 0.0)
  220. { 
  221. delete[] pnJs;
  222. return FALSE;
  223. }
  225. // 求解
  226.     pDataConst[n-1]=pDataConst[n-1]/d;
  227.     for (i=n-2;i>=0;i--)
  228.     { 
  229. t=0.0;
  230.         for (j=i+1;j<=n-1;j++)
  231. t=t+pDataCoef[i*n+j]*pDataConst[j];
  232.         pDataConst[i]=pDataConst[i]-t;
  233.     }
  234.     
  235. // 调整解的位置
  236. pnJs[n-1]=n-1;
  237.     for (k=n-1;k>=0;k--)
  238.     {
  239. if (pnJs[k]!=k)
  240.         { 
  241. t=pDataConst[k]; 
  242. pDataConst[k]=pDataConst[pnJs[k]]; 
  243. pDataConst[pnJs[k]]=t;
  244. }
  245. }
  247. // 清理内存
  248.     delete[] pnJs;
  250.     return TRUE;
  251. }
  253. //////////////////////////////////////////////////////////////////////
  254. // 全选主元高斯－约当消去法
  255. //
  256. // 参数：
  257. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组的解
  258. //
  259. // 返回值：BOOL 型，方程组求解是否成功
  260. //////////////////////////////////////////////////////////////////////
  261. BOOL CLEquations::GetRootsetGaussJordan(CMatrix& mtxResult)
  262. { 
  263. int *pnJs,l,k,i,j,nIs,p,q;
  264.     double d,t;
  266. // 方程组的属性，将常数矩阵赋给解矩阵
  267. mtxResult = m_mtxConst;
  268. double *pDataCoef = m_mtxCoef.GetData();
  269. double *pDataConst = mtxResult.GetData();
  270. int n = GetNumUnknowns();
  271. int m = m_mtxConst.GetNumColumns();
  273. // 临时缓冲区，存放变换的列数
  274.     pnJs = new int[n];
  276. // 消元
  277.     l=1;
  278.     for (k=0;k<=n-1;k++)
  279.     { 
  280. d=0.0;
  281.         for (i=k;i<=n-1;i++)
  282.         {
  283. for (j=k;j<=n-1;j++)
  284.             { 
  285. t=fabs(pDataCoef[i*n+j]);
  286. if (t>d) 
  287. { 
  288. d=t; 
  289. pnJs[k]=j; 
  290. nIs=i;
  291. }
  292.             }
  293. }
  295.         if (d+1.0==1.0) 
  296. l=0;
  297.         else
  298.         { 
  299. if (pnJs[k]!=k)
  300.             {
  301. for (i=0;i<=n-1;i++)
  302.                 { 
  303. p=i*n+k; 
  304. q=i*n+pnJs[k];
  305. t=pDataCoef[p]; 
  306. pDataCoef[p]=pDataCoef[q]; 
  307. pDataCoef[q]=t;
  308.                 }
  309. }
  311.             if (nIs!=k)
  312.             { 
  313. for (j=k;j<=n-1;j++)
  314.                 { 
  315. p=k*n+j; 
  316. q=nIs*n+j;
  317.                     t=pDataCoef[p]; 
  318. pDataCoef[p]=pDataCoef[q]; 
  319. pDataCoef[q]=t;
  320.                 }
  321.                 
  322. for (j=0;j<=m-1;j++)
  323.                 { 
  324. p=k*m+j; 
  325. q=nIs*m+j;
  326.                     t=pDataConst[p]; 
  327. pDataConst[p]=pDataConst[q]; 
  328. pDataConst[q]=t;
  329.                 }
  330.             }
  331.         }
  332.         
  333. // 求解失败
  334. if (l==0)
  335.         { 
  336. delete[] pnJs;
  337.             return FALSE;
  338.         }
  339.         
  340. d=pDataCoef[k*n+k];
  341.         for (j=k+1;j<=n-1;j++)
  342.         { 
  343. p=k*n+j; 
  344. pDataCoef[p]=pDataCoef[p]/d;
  345. }
  346.         
  347. for (j=0;j<=m-1;j++)
  348.         { 
  349. p=k*m+j; 
  350. pDataConst[p]=pDataConst[p]/d;
  351. }
  352.         
  353. for (j=k+1;j<=n-1;j++)
  354.         {
  355. for (i=0;i<=n-1;i++)
  356.             { 
  357. p=i*n+j;
  358. if (i!=k)
  359. pDataCoef[p]=pDataCoef[p]-pDataCoef[i*n+k]*pDataCoef[k*n+j];
  360.             }
  361. }
  363.         for (j=0;j<=m-1;j++)
  364.         {
  365. for (i=0;i<=n-1;i++)
  366. { 
  367. p=i*m+j;
  368. if (i!=k)
  369. pDataConst[p]=pDataConst[p]-pDataCoef[i*n+k]*pDataConst[k*m+j];
  370. }
  371. }
  372.     }
  373.     
  374. // 调整
  375. for (k=n-1;k>=0;k--)
  376.     {
  377. if (pnJs[k]!=k)
  378. {
  379. for (j=0;j<=m-1;j++)
  380. { 
  381. p=k*m+j; 
  382. q=pnJs[k]*m+j;
  383. t=pDataConst[p]; 
  384. pDataConst[p]=pDataConst[q]; 
  385. pDataConst[q]=t;
  386. }
  387. }
  388. }
  390. // 清理内存
  391.     delete[] pnJs;
  392.     
  393. return TRUE;
  394. }
  396. //////////////////////////////////////////////////////////////////////
  397. // 复系数方程组的全选主元高斯消去法
  398. //
  399. // 参数：
  400. // 1. const CMatrix& mtxCoefImag - 系数矩阵的虚部矩阵
  401. // 2. const CMatrix& mtxConstImag - 常数矩阵的虚部矩阵
  402. // 3. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵的实部矩阵
  403. // 4. CMatrix& mtxResultImag - CMatrix引用对象，返回方程组解矩阵的虚部矩阵
  404. //
  405. // 返回值：BOOL 型，方程组求解是否成功
  406. //////////////////////////////////////////////////////////////////////
  407. BOOL CLEquations::GetRootsetGauss(const CMatrix& mtxCoefImag, const CMatrix& mtxConstImag, CMatrix& mtxResult, CMatrix& mtxResultImag)
  408. { 
  409. int *pnJs,l,k,i,j,nIs,u,v;
  410.     double p,q,s,d;
  412. // 方程组的属性，将常数矩阵赋给解矩阵
  413. mtxResult = m_mtxConst;
  414. mtxResultImag = mtxConstImag;
  415. double *pDataCoef = m_mtxCoef.GetData();
  416. double *pDataConst = mtxResult.GetData();
  417. double *pDataCoefImag = mtxCoefImag.GetData();
  418. double *pDataConstImag = mtxResultImag.GetData();
  419. int n = GetNumUnknowns();
  420. int m = m_mtxConst.GetNumColumns();
  422. // 临时缓冲区，存放变换的列数
  423.     pnJs = new int[n];
  424.     
  425. // 消元
  426. for (k=0;k<=n-2;k++)
  427.     { 
  428. d=0.0;
  429.         for (i=k;i<=n-1;i++)
  430.         {
  431. for (j=k;j<=n-1;j++)
  432. { 
  433. u=i*n+j;
  434. p=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
  435. if (p>d) 
  436. {
  437. d=p;
  438. pnJs[k]=j;
  439. nIs=i;
  440. }
  441. }
  442.         }
  443.         
  444. // 求解失败
  445. if (d == 0.0)
  446. {
  447. delete[] pnJs;
  448. return FALSE;
  449.         }
  450.         
  451. if (nIs!=k)
  452.         { 
  453. for (j=k;j<=n-1;j++)
  454.             { 
  455. u=k*n+j; 
  456. v=nIs*n+j;
  457.                 p=pDataCoef[u]; 
  458. pDataCoef[u]=pDataCoef[v]; 
  459. pDataCoef[v]=p;
  460.                 p=pDataCoefImag[u]; 
  461. pDataCoefImag[u]=pDataCoefImag[v]; 
  462. pDataCoefImag[v]=p;
  463.             }
  464.             
  465. p=pDataConst[k]; 
  466. pDataConst[k]=pDataConst[nIs]; 
  467. pDataConst[nIs]=p;
  468.             p=pDataConstImag[k]; 
  469. pDataConstImag[k]=pDataConstImag[nIs]; 
  470. pDataConstImag[nIs]=p;
  471.         }
  472.         
  473. if (pnJs[k]!=k)
  474.         {
  475. for (i=0;i<=n-1;i++)
  476.             { 
  477. u=i*n+k; 
  478. v=i*n+pnJs[k];
  479. p=pDataCoef[u]; 
  480. pDataCoef[u]=pDataCoef[v]; 
  481. pDataCoef[v]=p;
  482. p=pDataCoefImag[u]; 
  483. pDataCoefImag[u]=pDataCoefImag[v]; 
  484. pDataCoefImag[v]=p;
  485.             }
  486. }
  488.         v=k*n+k;
  489.         for (j=k+1;j<=n-1;j++)
  490.         { 
  491. u=k*n+j;
  492.             p=pDataCoef[u]*pDataCoef[v]; 
  493. q=-pDataCoefImag[u]*pDataCoefImag[v];
  494.             s=(pDataCoef[v]-pDataCoefImag[v])*(pDataCoef[u]+pDataCoefImag[u]);
  495.             pDataCoef[u]=(p-q)/d; 
  496. pDataCoefImag[u]=(s-p-q)/d;
  497.         }
  498.         
  499. p=pDataConst[k]*pDataCoef[v]; 
  500. q=-pDataConstImag[k]*pDataCoefImag[v];
  501.         s=(pDataCoef[v]-pDataCoefImag[v])*(pDataConst[k]+pDataConstImag[k]);
  502.         pDataConst[k]=(p-q)/d; 
  503. pDataConstImag[k]=(s-p-q)/d;
  505.         for (i=k+1;i<=n-1;i++)
  506.         { 
  507. u=i*n+k;
  508.             for (j=k+1;j<=n-1;j++)
  509.             { 
  510. v=k*n+j; 
  511. l=i*n+j;
  512.                 p=pDataCoef[u]*pDataCoef[v]; 
  513. q=pDataCoefImag[u]*pDataCoefImag[v];
  514.                 s=(pDataCoef[u]+pDataCoefImag[u])*(pDataCoef[v]+pDataCoefImag[v]);
  515.                 pDataCoef[l]=pDataCoef[l]-p+q;
  516.                 pDataCoefImag[l]=pDataCoefImag[l]-s+p+q;
  517.             }
  518.             
  519. p=pDataCoef[u]*pDataConst[k]; 
  520. q=pDataCoefImag[u]*pDataConstImag[k];
  521.             s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[k]+pDataConstImag[k]);
  522.             pDataConst[i]=pDataConst[i]-p+q; 
  523. pDataConstImag[i]=pDataConstImag[i]-s+p+q;
  524.         }
  525.     }
  526.     
  527. u=(n-1)*n+n-1;
  528.     d=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
  530. // 求解失败
  531.     if (d == 0.0)
  532. {
  533. delete[] pnJs;
  534. return FALSE;
  535.     }
  537. // 求解
  538.     p=pDataCoef[u]*pDataConst[n-1]; q=-pDataCoefImag[u]*pDataConstImag[n-1];
  539.     s=(pDataCoef[u]-pDataCoefImag[u])*(pDataConst[n-1]+pDataConstImag[n-1]);
  540.     pDataConst[n-1]=(p-q)/d; pDataConstImag[n-1]=(s-p-q)/d;
  542.     for (i=n-2;i>=0;i--)
  543.     {
  544. for (j=i+1;j<=n-1;j++)
  545. { 
  546. u=i*n+j;
  547. p=pDataCoef[u]*pDataConst[j]; 
  548. q=pDataCoefImag[u]*pDataConstImag[j];
  549. s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[j]+pDataConstImag[j]);
  550. pDataConst[i]=pDataConst[i]-p+q;
  551. pDataConstImag[i]=pDataConstImag[i]-s+p+q;
  552. }
  553.     }
  555. // 调整位置
  556.     pnJs[n-1]=n-1;
  557.     for (k=n-1;k>=0;k--)
  558.     {
  559. if (pnJs[k]!=k)
  560.         { 
  561. p=pDataConst[k]; 
  562. pDataConst[k]=pDataConst[pnJs[k]]; 
  563. pDataConst[pnJs[k]]=p;
  564. p=pDataConstImag[k]; 
  565. pDataConstImag[k]=pDataConstImag[pnJs[k]]; 
  566. pDataConstImag[pnJs[k]]=p;
  567.         }
  568. }
  570. // 清理内存
  571. delete[] pnJs;
  573. return TRUE;
  574. }
  576. //////////////////////////////////////////////////////////////////////
  577. // 复系数方程组的全选主元高斯－约当消去法
  578. //
  579. // 参数：
  580. // 1. const CMatrix& mtxCoefImag - 系数矩阵的虚部矩阵
  581. // 2. const CMatrix& mtxConstImag - 常数矩阵的虚部矩阵
  582. // 3. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵的实部矩阵
  583. // 4. CMatrix& mtxResultImag - CMatrix引用对象，返回方程组解矩阵的虚部矩阵
  584. //
  585. // 返回值：BOOL 型，方程组求解是否成功
  586. //////////////////////////////////////////////////////////////////////
  587. BOOL CLEquations::GetRootsetGaussJordan(const CMatrix& mtxCoefImag, const CMatrix& mtxConstImag, CMatrix& mtxResult, CMatrix& mtxResultImag)
  588. {
  589. int *pnJs,l,k,i,j,nIs,u,v;
  590.     double p,q,s,d;
  592. // 方程组的属性，将常数矩阵赋给解矩阵
  593. mtxResult = m_mtxConst;
  594. mtxResultImag = mtxConstImag;
  595. double *pDataCoef = m_mtxCoef.GetData();
  596. double *pDataConst = mtxResult.GetData();
  597. double *pDataCoefImag = mtxCoefImag.GetData();
  598. double *pDataConstImag = mtxResultImag.GetData();
  599. int n = GetNumUnknowns();
  600. int m = m_mtxConst.GetNumColumns();
  602. // 临时缓冲区，存放变换的列数
  603.     pnJs = new int[n];
  605. // 消元
  606.     for (k=0;k<=n-1;k++)
  607.     { 
  608. d=0.0;
  609.         for (i=k;i<=n-1;i++)
  610.         {
  611. for (j=k;j<=n-1;j++)
  612. { 
  613. u=i*n+j;
  614. p=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
  615. if (p>d) 
  616. {
  617. d=p;
  618. pnJs[k]=j;
  619. nIs=i;
  620. }
  621. }
  622.         }
  623.         
  624. // 求解失败
  625. if (d == 0.0)
  626.         {
  627. delete[] pnJs;
  628.             return FALSE;
  629.         }
  630.         
  631. if (nIs!=k)
  632.         { 
  633. for (j=k;j<=n-1;j++)
  634.             { 
  635. u=k*n+j; 
  636. v=nIs*n+j;
  637.                 p=pDataCoef[u]; 
  638. pDataCoef[u]=pDataCoef[v]; 
  639. pDataCoef[v]=p;
  640.                 p=pDataCoefImag[u]; 
  641. pDataCoefImag[u]=pDataCoefImag[v]; 
  642. pDataCoefImag[v]=p;
  643.             }
  644.             
  645. for (j=0;j<=m-1;j++)
  646.             { 
  647. u=k*m+j; 
  648. v=nIs*m+j;
  649.                 p=pDataConst[u]; 
  650. pDataConst[u]=pDataConst[v]; 
  651. pDataConst[v]=p;
  652.                 p=pDataConstImag[u]; 
  653. pDataConstImag[u]=pDataConstImag[v]; 
  654. pDataConstImag[v]=p;
  655.             }
  656.         }
  657.         
  658. if (pnJs[k]!=k)
  659.         {
  660. for (i=0;i<=n-1;i++)
  661.             { 
  662. u=i*n+k; 
  663. v=i*n+pnJs[k];
  664. p=pDataCoef[u]; 
  665. pDataCoef[u]=pDataCoef[v]; 
  666. pDataCoef[v]=p;
  667. p=pDataCoefImag[u]; 
  668. pDataCoefImag[u]=pDataCoefImag[v]; 
  669. pDataCoefImag[v]=p;
  670.             }
  671. }
  673.         v=k*n+k;
  674.         for (j=k+1;j<=n-1;j++)
  675.         { 
  676. u=k*n+j;
  677.             p=pDataCoef[u]*pDataCoef[v]; 
  678. q=-pDataCoefImag[u]*pDataCoefImag[v];
  679.             s=(pDataCoef[v]-pDataCoefImag[v])*(pDataCoef[u]+pDataCoefImag[u]);
  680.             pDataCoef[u]=(p-q)/d; 
  681. pDataCoefImag[u]=(s-p-q)/d;
  682.         }
  683.         
  684. for (j=0;j<=m-1;j++)
  685.         { 
  686. u=k*m+j;
  687.             p=pDataConst[u]*pDataCoef[v]; 
  688. q=-pDataConstImag[u]*pDataCoefImag[v];
  689.             s=(pDataCoef[v]-pDataCoefImag[v])*(pDataConst[u]+pDataConstImag[u]);
  690.             pDataConst[u]=(p-q)/d; 
  691. pDataConstImag[u]=(s-p-q)/d;
  692.         }
  693.         
  694. for (i=0;i<=n-1;i++)
  695.         {
  696. if (i!=k)
  697.             { 
  698. u=i*n+k;
  699. for (j=k+1;j<=n-1;j++)
  700.                 { 
  701. v=k*n+j; 
  702. l=i*n+j;
  703. p=pDataCoef[u]*pDataCoef[v]; 
  704. q=pDataCoefImag[u]*pDataCoefImag[v];
  705. s=(pDataCoef[u]+pDataCoefImag[u])*(pDataCoef[v]+pDataCoefImag[v]);
  706. pDataCoef[l]=pDataCoef[l]-p+q;
  707. pDataCoefImag[l]=pDataCoefImag[l]-s+p+q;
  708.                 }
  709.             
  710. for (j=0;j<=m-1;j++)
  711. { 
  712. l=i*m+j; 
  713. v=k*m+j;
  714. p=pDataCoef[u]*pDataConst[v]; q=pDataCoefImag[u]*pDataConstImag[v];
  715. s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[v]+pDataConstImag[v]);
  716. pDataConst[l]=pDataConst[l]-p+q; 
  717. pDataConstImag[l]=pDataConstImag[l]-s+p+q;
  718. }
  719.             }
  720. }
  721.     }
  723. // 求解调整
  724. for (k=n-1;k>=0;k--)
  725.     {
  726. if (pnJs[k]!=k)
  727. {
  728. for (j=0;j<=m-1;j++)
  729. { 
  730. u=k*m+j;
  731. v=pnJs[k]*m+j;
  732. p=pDataConst[u]; 
  733. pDataConst[u]=pDataConst[v]; 
  734. pDataConst[v]=p;
  735. p=pDataConstImag[u]; 
  736. pDataConstImag[u]=pDataConstImag[v]; 
  737. pDataConstImag[v]=p;
  738. }
  739. }
  740. }
  742. // 清理内存
  743. delete[] pnJs;
  745. return TRUE;
  746. }
  748. //////////////////////////////////////////////////////////////////////
  749. // 求解三对角线方程组的追赶法
  750. //
  751. // 参数：
  752. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  753. //
  754. // 返回值：BOOL 型，方程组求解是否成功
  755. //////////////////////////////////////////////////////////////////////
  756. BOOL CLEquations::GetRootsetTriDiagonal(CMatrix& mtxResult)
  757. { 
  758. int k,j;
  759.     double s;
  760.     
  761. // 将常数矩阵赋给解矩阵
  762. mtxResult = m_mtxConst;
  763. double *pDataConst = mtxResult.GetData();
  765. int n = GetNumUnknowns();
  766. ASSERT(m_mtxCoef.GetNumRows() == n);
  768. // 为系数矩阵对角线数组分配内存
  769. double* pDiagData = new double[3*n-2];
  771. // 构造系数矩阵对角线元素数组
  772. k = j = 0;
  773. if (k == 0)
  774. {
  775. pDiagData[j++] = m_mtxCoef.GetElement(k, k);
  776. pDiagData[j++] = m_mtxCoef.GetElement(k, k+1);
  777. }
  778.     for (k=1; k<n-1; ++k)
  779. {
  780. pDiagData[j++] = m_mtxCoef.GetElement(k, k-1);
  781. pDiagData[j++] = m_mtxCoef.GetElement(k, k);
  782. pDiagData[j++] = m_mtxCoef.GetElement(k, k+1);
  783. }
  784. if (k == n-1)
  785. {
  786. pDiagData[j++] = m_mtxCoef.GetElement(k, k-1);
  787. pDiagData[j++] = m_mtxCoef.GetElement(k, k);
  788. }
  790. // 求解
  791. for (k=0;k<=n-2;k++)
  792.     { 
  793. j=3*k; 
  794. s=pDiagData[j];
  796. // 求解失败
  797.         if (fabs(s)+1.0==1.0)
  798. {
  799. delete[] pDiagData;
  800. return FALSE;
  801. }
  803.         pDiagData[j+1]=pDiagData[j+1]/s;
  804.         pDataConst[k]=pDataConst[k]/s;
  805.         pDiagData[j+3]=pDiagData[j+3]-pDiagData[j+2]*pDiagData[j+1];
  806.         pDataConst[k+1]=pDataConst[k+1]-pDiagData[j+2]*pDataConst[k];
  807.     }
  808.     
  809. s=pDiagData[3*n-3];
  810.     if (s == 0.0)
  811. {
  812. delete[] pDiagData;
  813. return FALSE;
  814. }
  815.     
  816. // 调整
  817. pDataConst[n-1]=pDataConst[n-1]/s;
  818.     for (k=n-2;k>=0;k--)
  819. pDataConst[k]=pDataConst[k]-pDiagData[3*k+1]*pDataConst[k+1];
  820.     
  821. // 释放内存
  822. delete[] pDiagData;
  824. return TRUE;
  825. }
  827. //////////////////////////////////////////////////////////////////////
  828. // 一般带型方程组的求解
  829. //
  830. // 参数：
  831. // 1. int nBandWidth - 带宽
  832. // 2. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  833. //
  834. // 返回值：BOOL 型，方程组求解是否成功
  835. //////////////////////////////////////////////////////////////////////
  836. BOOL CLEquations::GetRootsetBand(int nBandWidth, CMatrix& mtxResult)
  837. { 
  838. int ls,k,i,j,is,u,v;
  839.     double p,t;
  840.     
  841. // 带宽必须为奇数
  842. if ((nBandWidth-1)%2 != 0)
  843. return FALSE;
  845. // 将常数矩阵赋给解矩阵
  846. mtxResult = m_mtxConst;
  847. double *pDataConst = mtxResult.GetData();
  849. // 方程组属性
  850. int m = m_mtxConst.GetNumColumns();
  851. int n = GetNumUnknowns();
  852. ASSERT(m_mtxCoef.GetNumRows() == n);
  854. // 带宽数组：带型矩阵的有效数据
  855. double* pBandData = new double[nBandWidth*n];
  857. // 半带宽
  858.     ls = (nBandWidth-1)/2;
  860. // 构造带宽数组
  861.     for (i=0; i<n; ++i)
  862. {
  863. j = 0;
  864. for (k=max(0, i-ls); k<max(0, i-ls)+nBandWidth; ++k)
  865. {
  866. if (k < n)
  867. pBandData[i*nBandWidth+j++] = m_mtxCoef.GetElement(i, k);
  868. else
  869. pBandData[i*nBandWidth+j++] = 0;
  870. }
  871. }
  873. // 求解
  874.     for (k=0;k<=n-2;k++)
  875.     { 
  876. p=0.0;
  877.         for (i=k;i<=ls;i++)
  878.         { 
  879. t=fabs(pBandData[i*nBandWidth]);
  880.             if (t>p) 
  881. {
  882. p=t; 
  883. is=i;
  884. }
  885.         }
  886.         
  887. if (p == 0.0)
  888. {
  889. delete[] pBandData;
  890. return FALSE;
  891. }
  893.         for (j=0;j<=m-1;j++)
  894.         { 
  895. u=k*m+j; 
  896. v=is*m+j;
  897.             t=pDataConst[u]; 
  898. pDataConst[u]=pDataConst[v]; 
  899. pDataConst[v]=t;
  900.         }
  901.         
  902. for (j=0;j<=nBandWidth-1;j++)
  903.         { 
  904. u=k*nBandWidth+j; 
  905. v=is*nBandWidth+j;
  906.             t=pBandData[u]; 
  907. pBandData[u]=pBandData[v]; 
  908. pBandData[v]=t;
  909.         }
  910.         
  911. for (j=0;j<=m-1;j++)
  912.         { 
  913. u=k*m+j; 
  914. pDataConst[u]=pDataConst[u]/pBandData[k*nBandWidth];
  915. }
  916.         
  917. for (j=1;j<=nBandWidth-1;j++)
  918.         { 
  919. u=k*nBandWidth+j; 
  920. pBandData[u]=pBandData[u]/pBandData[k*nBandWidth];
  921. }
  922.         
  923. for (i=k+1;i<=ls;i++)
  924.         { 
  925. t=pBandData[i*nBandWidth];
  926.             for (j=0;j<=m-1;j++)
  927.             { 
  928. u=i*m+j; 
  929. v=k*m+j;
  930.                 pDataConst[u]=pDataConst[u]-t*pDataConst[v];
  931.             }
  932.             
  933. for (j=1;j<=nBandWidth-1;j++)
  934.             { 
  935. u=i*nBandWidth+j; 
  936. v=k*nBandWidth+j;
  937.                 pBandData[u-1]=pBandData[u]-t*pBandData[v];
  938.             }
  939.             
  940. u=i*nBandWidth+nBandWidth-1; pBandData[u]=0.0;
  941.         }
  942.         
  943. if (ls!=(n-1)) 
  944. ls=ls+1;
  945.     }
  946.     
  947. p=pBandData[(n-1)*nBandWidth];
  948.     if (p == 0.0)
  949. {
  950. delete[] pBandData;
  951. return FALSE;
  952. }
  954.     for (j=0;j<=m-1;j++)
  955.     { 
  956. u=(n-1)*m+j; 
  957. pDataConst[u]=pDataConst[u]/p;
  958. }
  959.     
  960. ls=1;
  961.     for (i=n-2;i>=0;i--)
  962.     { 
  963. for (k=0;k<=m-1;k++)
  964.         { 
  965. u=i*m+k;
  966.             for (j=1;j<=ls;j++)
  967.             { 
  968. v=i*nBandWidth+j; 
  969. is=(i+j)*m+k;
  970.                 pDataConst[u]=pDataConst[u]-pBandData[v]*pDataConst[is];
  971.             }
  972.         }
  973.         
  974. if (ls!=(nBandWidth-1)) 
  975. ls=ls+1;
  976.     }
  977.     
  978. // 释放内存
  979. delete[] pBandData;
  981. return TRUE;
  982. }
  984. //////////////////////////////////////////////////////////////////////
  985. // 求解对称方程组的分解法
  986. //
  987. // 参数：
  988. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  989. //
  990. // 返回值：BOOL 型，方程组求解是否成功
  991. //////////////////////////////////////////////////////////////////////
  992. BOOL CLEquations::GetRootsetDjn(CMatrix& mtxResult)
  993. { 
  994. int i,j,l,k,u,v,w,k1,k2,k3;
  995.     double p;
  997. // 方程组属性，将常数矩阵赋给解矩阵
  998. CMatrix mtxCoef = m_mtxCoef;
  999. mtxResult = m_mtxConst;
  1000. int n = mtxCoef.GetNumColumns();
  1001. int m = mtxResult.GetNumColumns();
  1002. double* pDataCoef = mtxCoef.GetData();
  1003. double* pDataConst = mtxResult.GetData();
  1005. // 非对称系数矩阵，不能用本方法求解
  1006.     if (pDataCoef[0] == 0.0)
  1007. return FALSE;
  1009.     for (i=1; i<=n-1; i++)
  1010.     { 
  1011. u=i*n; 
  1012. pDataCoef[u]=pDataCoef[u]/pDataCoef[0];
  1013. }
  1014.     
  1015. for (i=1; i<=n-2; i++)
  1016.     { 
  1017. u=i*n+i;
  1018.         for (j=1; j<=i; j++)
  1019.         { 
  1020. v=i*n+j-1; 
  1021. l=(j-1)*n+j-1;
  1022.             pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v]*pDataCoef[l];
  1023.         }
  1024.         
  1025. p=pDataCoef[u];
  1026.         if (p == 0.0)
  1027. return FALSE;
  1029.         for (k=i+1; k<=n-1; k++)
  1030.         { 
  1031. u=k*n+i;
  1032.             for (j=1; j<=i; j++)
  1033.             { 
  1034. v=k*n+j-1; 
  1035. l=i*n+j-1; 
  1036. w=(j-1)*n+j-1;
  1037. pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[l]*pDataCoef[w];
  1038.             }
  1039.             
  1040. pDataCoef[u]=pDataCoef[u]/p;
  1041.         }
  1042.     }
  1043.     
  1044. u=n*n-1;
  1045.     for (j=1; j<=n-1; j++)
  1046.     { 
  1047. v=(n-1)*n+j-1; 
  1048. w=(j-1)*n+j-1;
  1049.         pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v]*pDataCoef[w];
  1050.     }
  1051.     
  1052. p=pDataCoef[u];
  1053.     if (p == 0.0)
  1054. return FALSE;
  1056.     for (j=0; j<=m-1; j++)
  1057. {
  1058.     for (i=1; i<=n-1; i++)
  1059. { 
  1060. u=i*m+j;
  1061. for (k=1; k<=i; k++)
  1062. { 
  1063. v=i*n+k-1; 
  1064. w=(k-1)*m+j;
  1065. pDataConst[u]=pDataConst[u]-pDataCoef[v]*pDataConst[w];
  1066. }
  1067. }
  1068. }
  1070.     for (i=1; i<=n-1; i++)
  1071.     { 
  1072. u=(i-1)*n+i-1;
  1073.         for (j=i; j<=n-1; j++)
  1074.         { 
  1075. v=(i-1)*n+j; 
  1076. w=j*n+i-1;
  1077.             pDataCoef[v]=pDataCoef[u]*pDataCoef[w];
  1078.         }
  1079.     }
  1080.     
  1081. for (j=0; j<=m-1; j++)
  1082.     { 
  1083. u=(n-1)*m+j;
  1084.         pDataConst[u]=pDataConst[u]/p;
  1085.         for (k=1; k<=n-1; k++)
  1086.         { 
  1087. k1=n-k; 
  1088. k3=k1-1; 
  1089. u=k3*m+j;
  1090.             for (k2=k1; k2<=n-1; k2++)
  1091.             { 
  1092. v=k3*n+k2; 
  1093. w=k2*m+j;
  1094.                 pDataConst[u]=pDataConst[u]-pDataCoef[v]*pDataConst[w];
  1095.             }
  1096.             
  1097. pDataConst[u]=pDataConst[u]/pDataCoef[k3*n+k3];
  1098.         }
  1099.     }
  1100.     
  1101. return TRUE;
  1102. }
  1104. //////////////////////////////////////////////////////////////////////
  1105. // 求解对称正定方程组的平方根法
  1106. //
  1107. // 参数：
  1108. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1109. //
  1110. // 返回值：BOOL 型，方程组求解是否成功
  1111. //////////////////////////////////////////////////////////////////////
  1112. BOOL CLEquations::GetRootsetCholesky(CMatrix& mtxResult)
  1113. { 
  1114. int i,j,k,u,v;
  1115.     
  1116. // 方程组属性，将常数矩阵赋给解矩阵
  1117. CMatrix mtxCoef = m_mtxCoef;
  1118. mtxResult = m_mtxConst;
  1119. int n = mtxCoef.GetNumColumns();
  1120. int m = mtxResult.GetNumColumns();
  1121. double* pDataCoef = mtxCoef.GetData();
  1122. double* pDataConst = mtxResult.GetData();
  1123.     
  1124. // 非对称正定系数矩阵，不能用本方法求解
  1125. if (pDataCoef[0] <= 0.0)
  1126. return FALSE;
  1128. pDataCoef[0]=sqrt(pDataCoef[0]);
  1129.     for (j=1; j<=n-1; j++) 
  1130. pDataCoef[j]=pDataCoef[j]/pDataCoef[0];
  1131.     
  1132. for (i=1; i<=n-1; i++)
  1133.     { 
  1134. u=i*n+i;
  1135.         for (j=1; j<=i; j++)
  1136.         { 
  1137. v=(j-1)*n+i;
  1138.             pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v];
  1139.         }
  1140.         
  1141. if (pDataCoef[u] <= 0.0)
  1142. return FALSE;
  1144.         pDataCoef[u]=sqrt(pDataCoef[u]);
  1145.         if (i!=(n-1))
  1146.         { 
  1147. for (j=i+1; j<=n-1; j++)
  1148.             { 
  1149. v=i*n+j;
  1150.                 for (k=1; k<=i; k++)
  1151. pDataCoef[v]=pDataCoef[v]-pDataCoef[(k-1)*n+i]*pDataCoef[(k-1)*n+j];
  1152.                 pDataCoef[v]=pDataCoef[v]/pDataCoef[u];
  1153.             }
  1154.         }
  1155.     }
  1156.     
  1157. for (j=0; j<=m-1; j++)
  1158.     { 
  1159. pDataConst[j]=pDataConst[j]/pDataCoef[0];
  1160.         for (i=1; i<=n-1; i++)
  1161.         { 
  1162. u=i*n+i; 
  1163. v=i*m+j;
  1164.             for (k=1; k<=i; k++)
  1165. pDataConst[v]=pDataConst[v]-pDataCoef[(k-1)*n+i]*pDataConst[(k-1)*m+j];
  1166.             pDataConst[v]=pDataConst[v]/pDataCoef[u];
  1167.         }
  1168.     }
  1169.     
  1170. for (j=0; j<=m-1; j++)
  1171.     { 
  1172. u=(n-1)*m+j;
  1173.         pDataConst[u]=pDataConst[u]/pDataCoef[n*n-1];
  1174.         for (k=n-1; k>=1; k--)
  1175.         { 
  1176. u=(k-1)*m+j;
  1177.             for (i=k; i<=n-1; i++)
  1178.             { 
  1179. v=(k-1)*n+i;
  1180.                 pDataConst[u]=pDataConst[u]-pDataCoef[v]*pDataConst[i*m+j];
  1181.             }
  1182.             
  1183. v=(k-1)*n+k-1;
  1184.             pDataConst[u]=pDataConst[u]/pDataCoef[v];
  1185.         }
  1186.     }
  1187.     
  1188. return TRUE;
  1189. }
  1191. //////////////////////////////////////////////////////////////////////
  1192. // 求解大型稀疏方程组的全选主元高斯－约去消去法
  1193. //
  1194. // 参数：
  1195. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1196. //
  1197. // 返回值：BOOL 型，方程组求解是否成功
  1198. //////////////////////////////////////////////////////////////////////
  1199. BOOL CLEquations::GetRootsetGgje(CMatrix& mtxResult)
  1200. { 
  1201. int *pnJs,i,j,k,nIs,u,v;
  1202.     double d,t;
  1203.     
  1204. // 方程组属性，将常数矩阵赋给解矩阵
  1205. CMatrix mtxCoef = m_mtxCoef;
  1206. mtxResult = m_mtxConst;
  1207. int n = mtxCoef.GetNumColumns();
  1208. double* pDataCoef = mtxCoef.GetData();
  1209. double* pDataConst = mtxResult.GetData();
  1211. // 临时缓冲区，存放变换的列数
  1212. pnJs = new int[n];
  1214. // 消元
  1215.     for (k=0; k<=n-1; k++)
  1216.     { 
  1217. d=0.0;
  1218.         for (i=k; i<=n-1; i++)
  1219. {
  1220. for (j=k; j<=n-1; j++)
  1221. { 
  1222. t=fabs(pDataCoef[i*n+j]);
  1223. if (t>d) 
  1224. {
  1225. d=t; 
  1226. pnJs[k]=j; 
  1227. nIs=i;
  1228. }
  1229. }
  1230. }
  1232.         if (d == 0.0)
  1233. {
  1234. delete[] pnJs;
  1235. return FALSE;
  1236. }
  1238.         if (nIs!=k)
  1239.         { 
  1240. for (j=k; j<=n-1; j++)
  1241.             { 
  1242. u=k*n+j; 
  1243. v=nIs*n+j;
  1244.                 t=pDataCoef[u]; 
  1245. pDataCoef[u]=pDataCoef[v]; 
  1246. pDataCoef[v]=t;
  1247.             }
  1248.             
  1249. t=pDataConst[k]; 
  1250. pDataConst[k]=pDataConst[nIs]; 
  1251. pDataConst[nIs]=t;
  1252.         }
  1253.         
  1254. if (pnJs[k]!=k)
  1255.         {
  1256. for (i=0; i<=n-1; i++)
  1257.             { 
  1258. u=i*n+k; 
  1259. v=i*n+pnJs[k];
  1260. t=pDataCoef[u]; 
  1261. pDataCoef[u]=pDataCoef[v]; 
  1262. pDataCoef[v]=t;
  1263.             }
  1264. }
  1266.         t=pDataCoef[k*n+k];
  1267.         for (j=k+1; j<=n-1; j++)
  1268.         { 
  1269. u=k*n+j;
  1270.             if (pDataCoef[u]!=0.0) 
  1271. pDataCoef[u]=pDataCoef[u]/t;
  1272.         }
  1273.         
  1274. pDataConst[k]=pDataConst[k]/t;
  1275.         for (j=k+1; j<=n-1; j++)
  1276.         { 
  1277. u=k*n+j;
  1278.             if (pDataCoef[u]!=0.0)
  1279.             { 
  1280. for (i=0; i<=n-1; i++)
  1281.                 { 
  1282. v=i*n+k;
  1283.                     if ((i!=k)&&(pDataCoef[v]!=0.0))
  1284.                     { 
  1285. nIs=i*n+j;
  1286.                         pDataCoef[nIs]=pDataCoef[nIs]-pDataCoef[v]*pDataCoef[u];
  1287.                     }
  1288.                 }
  1289.             }
  1290.         }
  1291.         
  1292. for (i=0; i<=n-1; i++)
  1293.         { 
  1294. u=i*n+k;
  1295.             if ((i!=k)&&(pDataCoef[u]!=0.0))
  1296. pDataConst[i]=pDataConst[i]-pDataCoef[u]*pDataConst[k];
  1297.         }
  1298.     }
  1299.     
  1300. // 调整
  1301. for (k=n-1; k>=0; k--)
  1302.     {
  1303. if (k!=pnJs[k])
  1304.         { 
  1305. t=pDataConst[k]; 
  1306. pDataConst[k]=pDataConst[pnJs[k]]; 
  1307. pDataConst[pnJs[k]]=t;
  1308. }
  1309. }
  1311. // 释放内存
  1312.     delete[] pnJs;
  1313.     
  1314. return TRUE;
  1315. }
  1317. //////////////////////////////////////////////////////////////////////
  1318. // 求解托伯利兹方程组的列文逊方法
  1319. //
  1320. // 参数：
  1321. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1322. //
  1323. // 返回值：BOOL 型，方程组求解是否成功
  1324. //////////////////////////////////////////////////////////////////////
  1325. BOOL CLEquations::GetRootsetTlvs(CMatrix& mtxResult)
  1326. { 
  1327. int i,j,k;
  1328.     double a,beta,q,c,h,*y,*s;
  1330. // 未知数个数
  1331.     int n = m_mtxCoef.GetNumColumns();
  1333. // 初始化解解向量
  1334. mtxResult.Init(n, 1);
  1335. double* x = mtxResult.GetData();
  1337. // 常数数组
  1338. double* pDataConst = m_mtxConst.GetData();
  1340. // 建立T数组
  1341. double* t = new double[n];
  1343. // 构造T数组
  1344.     for (i=0; i<n; ++i)
  1345. t[i] = m_mtxCoef.GetElement(0, i);
  1347. // 临时数组
  1348. s = new double[n];
  1349.     y = new double[n];
  1351. // 非托伯利兹方程组，不能用本方法求解
  1352.     a=t[0];
  1353.     if (a == 0.0)
  1354.     { 
  1355. delete[] s;
  1356. delete[] y;
  1357. delete[] t;
  1358. return FALSE;
  1359. }
  1361. // 列文逊方法求解
  1362.     y[0]=1.0; 
  1363. x[0]=pDataConst[0]/a;
  1364.     for (k=1; k<=n-1; k++)
  1365.     { 
  1366. beta=0.0; 
  1367. q=0.0;
  1368.         for (j=0; j<=k-1; j++)
  1369.         { 
  1370. beta=beta+y[j]*t[j+1];
  1371.             q=q+x[j]*t[k-j];
  1372.         }
  1373.         
  1374. if (a == 0.0)
  1375. { 
  1376. delete[] s;
  1377. delete[] y;
  1378. delete[] t;
  1379. return FALSE;
  1380. }
  1382.         c=-beta/a; 
  1383. s[0]=c*y[k-1]; 
  1384. y[k]=y[k-1];
  1385.         if (k!=1)
  1386.         {
  1387. for (i=1; i<=k-1; i++)
  1388. s[i]=y[i-1]+c*y[k-i-1];
  1389. }
  1391.         a=a+c*beta;
  1392. if (a == 0.0)
  1393. { 
  1394. delete[] s;
  1395. delete[] y;
  1396. delete[] t;
  1397. return FALSE;
  1398. }
  1400.         h=(pDataConst[k]-q)/a;
  1401.         for (i=0; i<=k-1; i++)
  1402.         { 
  1403. x[i]=x[i]+h*s[i]; 
  1404. y[i]=s[i];
  1405. }
  1406.         
  1407. x[k]=h*y[k];
  1408.     }
  1409.     
  1410. // 释放内存
  1411. delete[] s;
  1412. delete[] y;
  1413. delete[] t;
  1415.     return TRUE;
  1416. }
  1418. //////////////////////////////////////////////////////////////////////
  1419. // 高斯－赛德尔迭代法
  1420. //
  1421. // 参数：
  1422. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1423. // 2. double eps - 控制精度，默认值为0.000001
  1424. //
  1425. // 返回值：BOOL 型，方程组求解是否成功
  1426. //////////////////////////////////////////////////////////////////////
  1427. BOOL CLEquations::GetRootsetGaussSeidel(CMatrix& mtxResult, double eps /*= 0.000001*/)
  1428. { 
  1429. int i,j,u,v;
  1430.     double p,t,s,q;
  1432. // 未知数个数
  1433.     int n = m_mtxCoef.GetNumColumns();
  1435. // 初始化解向量
  1436. mtxResult.Init(n, 1);
  1437. double* x = mtxResult.GetData();
  1439. // 系数与常数
  1440. double* pDataCoef = m_mtxCoef.GetData();
  1441. double* pDataConst = m_mtxConst.GetData();
  1442.     
  1443. // 求解
  1444. for (i=0; i<=n-1; i++)
  1445.     { 
  1446. u=i*n+i; 
  1447. p=0.0; 
  1448. x[i]=0.0;
  1449.         for (j=0; j<=n-1; j++)
  1450.         {
  1451. if (i!=j)
  1452.             { 
  1453. v=i*n+j; 
  1454. p=p+fabs(pDataCoef[v]);
  1455. }
  1456. }
  1458.         if (p>=fabs(pDataCoef[u]))
  1459. return FALSE;
  1460.     }
  1462. // 精度控制
  1463.     p=eps+1.0;
  1464.     while (p>=eps)
  1465.     { 
  1466. p=0.0;
  1467.         for (i=0; i<=n-1; i++)
  1468.         { 
  1469. t=x[i]; 
  1470. s=0.0;
  1471.             for (j=0; j<=n-1; j++)
  1472. if (j!=i) 
  1473. s=s+pDataCoef[i*n+j]*x[j];
  1475.             x[i]=(pDataConst[i]-s)/pDataCoef[i*n+i];
  1476.             q=fabs(x[i]-t)/(1.0+fabs(x[i]));
  1477.             if (q>p) 
  1478. p=q;
  1479.         }
  1480.     }
  1481.     
  1482. return TRUE;
  1483. }
  1485. //////////////////////////////////////////////////////////////////////
  1486. // 求解对称正定方程组的共轭梯度法
  1487. //
  1488. // 参数：
  1489. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1490. // 2. double eps - 控制精度，默认值为0.000001
  1491. //
  1492. // 返回值：BOOL 型，方程组求解是否成功
  1493. //////////////////////////////////////////////////////////////////////
  1494. void CLEquations::GetRootsetGrad(CMatrix& mtxResult, double eps /*= 0.000001*/)
  1495. { 
  1496. int i,k;
  1497.     double *p,*r,*s,*q,alpha,beta,d,e;
  1499. // 未知数个数
  1500.     int n = GetNumUnknowns();
  1502. // 初始化解向量
  1503. mtxResult.Init(n, 1);
  1504. double* x = mtxResult.GetData();
  1506. // 构造临时矩阵
  1507. CMatrix mtxP(n, 1);
  1508.     p = mtxP.GetData();
  1510. double* pDataCoef = m_mtxCoef.GetData();
  1511. double* pDataConst = m_mtxConst.GetData();
  1513.     r = new double[n];
  1515.     for (i=0; i<=n-1; i++)
  1516.     { 
  1517. x[i]=0.0; 
  1518. p[i]=pDataConst[i]; 
  1519. r[i]=pDataConst[i]; 
  1520. }
  1521.     
  1522. i=0;
  1523.     while (i<=n-1)
  1524.     { 
  1525. CMatrix mtxS = m_mtxCoef * mtxP;
  1526.     s = mtxS.GetData();
  1527.         
  1528. d=0.0; 
  1529. e=0.0;
  1530.         for (k=0; k<=n-1; k++)
  1531.         { 
  1532. d=d+p[k]*pDataConst[k]; 
  1533. e=e+p[k]*s[k]; 
  1534. }
  1535.         
  1536. alpha=d/e;
  1537.         for (k=0; k<=n-1; k++)
  1538. x[k]=x[k]+alpha*p[k];
  1539.         
  1540. CMatrix mtxQ = m_mtxCoef * mtxResult;
  1541.     q = mtxQ.GetData();
  1542.         
  1543. d=0.0;
  1544.         for (k=0; k<=n-1; k++)
  1545.         { 
  1546. r[k]=pDataConst[k]-q[k]; 
  1547. d=d+r[k]*s[k]; 
  1548. }
  1549.         
  1550. beta=d/e; d=0.0;
  1551.         for (k=0; k<=n-1; k++) 
  1552. d=d+r[k]*r[k];
  1553.         
  1554. // 满足精度，求解结束
  1555. d=sqrt(d);
  1556.         if (d<eps)
  1557. break;
  1559.         for (k=0; k<=n-1; k++)
  1560. p[k]=r[k]-beta*p[k];
  1561.         
  1562. i=i+1;
  1563.     }
  1565. delete[] r;
  1566. }
  1568. //////////////////////////////////////////////////////////////////////
  1569. // 求解线性最小二乘问题的豪斯荷尔德变换法
  1570. //
  1571. // 参数：
  1572. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1573. // 2. CMatrix& mtxQ - CMatrix引用对象，返回豪斯荷尔德变换的Q矩阵
  1574. // 3. CMatrix& mtxR - CMatrix引用对象，返回豪斯荷尔德变换的R矩阵
  1575. //
  1576. // 返回值：BOOL 型，方程组求解是否成功
  1577. //////////////////////////////////////////////////////////////////////
  1578. BOOL CLEquations::GetRootsetMqr(CMatrix& mtxResult, CMatrix& mtxQ, CMatrix& mtxR)
  1579. { 
  1580. int i,j;
  1581.     double d;
  1583. // 方程组的方程数和未知数个数
  1584.     int m = m_mtxCoef.GetNumRows();
  1585.     int n = m_mtxCoef.GetNumColumns();
  1586. // 奇异方程组
  1587. if (m < n)
  1588. return FALSE;
  1590. // 将解向量初始化为常数向量
  1591. mtxResult = m_mtxConst;
  1592. double* pDataConst = mtxResult.GetData();
  1594. // 构造临时矩阵，用于QR分解
  1595. mtxR = m_mtxCoef;
  1596. double* pDataCoef = mtxR.GetData();
  1598. // QR分解
  1599.     if (! mtxR.SplitQR(mtxQ))
  1600. return FALSE;
  1602. // 临时缓冲区
  1603. double* c = new double[n];
  1604. double* q = mtxQ.GetData();
  1606. // 求解
  1607.     for (i=0; i<=n-1; i++)
  1608.     { 
  1609. d=0.0;
  1610.         for (j=0; j<=m-1; j++)
  1611. d=d+q[j*m+i]*pDataConst[j];
  1612.     
  1613. c[i]=d;
  1614.     }
  1615.     
  1616. pDataConst[n-1]=c[n-1]/pDataCoef[n*n-1];
  1617.     for (i=n-2; i>=0; i--)
  1618.     { 
  1619. d=0.0;
  1620.         for (j=i+1; j<=n-1; j++)
  1621. d=d+pDataCoef[i*n+j]*pDataConst[j];
  1622.         
  1623. pDataConst[i]=(c[i]-d)/pDataCoef[i*n+i];
  1624.     }
  1625.     
  1626. // 释放内存
  1627. delete[] c;
  1629. return TRUE;
  1630. }
  1632. //////////////////////////////////////////////////////////////////////
  1633. // 求解线性最小二乘问题的广义逆法
  1634. //
  1635. // 参数：
  1636. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1637. // 2. CMatrix& mtxAP - CMatrix引用对象，返回系数矩阵的广义逆矩阵
  1638. // 3. CMatrix& mtxU - CMatrix引用对象，返回U矩阵
  1639. // 4. CMatrix& mtxV - CMatrix引用对象，返回V矩阵
  1640. // 5. double eps - 控制精度，默认值为0.000001
  1641. //
  1642. // 返回值：BOOL 型，方程组求解是否成功
  1643. //////////////////////////////////////////////////////////////////////
  1644. BOOL CLEquations::GetRootsetGinv(CMatrix& mtxResult, CMatrix& mtxAP, CMatrix& mtxU, CMatrix& mtxV, double eps /*= 0.000001*/)
  1645. { 
  1646. int i,j;
  1647.     
  1648. // 方程个数和未知数个数
  1649.     int m = m_mtxCoef.GetNumRows();
  1650.     int n = m_mtxCoef.GetNumColumns();
  1652. // 初始化解向量
  1653. mtxResult.Init(n, 1);
  1655. double* pDataConst = m_mtxConst.GetData();
  1656. double* x = mtxResult.GetData();
  1658. // 临时矩阵
  1659. CMatrix mtxA = m_mtxCoef;
  1661. // 求广义逆矩阵
  1662.     if (! mtxA.GInvertUV(mtxAP, mtxU, mtxV, eps))
  1663. return FALSE;
  1665. double* pAPData = mtxAP.GetData();
  1667. // 求解
  1668.     for (i=0; i<=n-1; i++)
  1669.     { 
  1670. x[i]=0.0;
  1671.         for (j=0; j<=m-1; j++)
  1672. x[i]=x[i]+pAPData[i*m+j]*pDataConst[j];
  1673.     }
  1674.     
  1675. return TRUE;
  1676. }
  1678. //////////////////////////////////////////////////////////////////////
  1679. // 病态方程组的求解
  1680. //
  1681. // 参数：
  1682. // 1. CMatrix& mtxResult - CMatrix引用对象，返回方程组解矩阵
  1683. // 2. int nMaxIt - 叠加次数，默认值为60
  1684. // 3. double eps - 控制精度，默认值为0.000001
  1685. //
  1686. // 返回值：BOOL 型，方程组求解是否成功
  1687. //////////////////////////////////////////////////////////////////////
  1688. BOOL CLEquations::GetRootsetMorbid(CMatrix& mtxResult, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  1689. { 
  1690. int i, k;
  1691.     double q, qq;
  1692.     
  1693. // 方程的阶数
  1694.     int n = GetNumUnknowns();
  1696. // 设定迭代次数, 缺省为60
  1697. i = nMaxIt;
  1698.     
  1699. // 用全选主元高斯消元法求解
  1700. CLEquations leqs(m_mtxCoef, m_mtxConst);
  1701. if (! leqs.GetRootsetGauss(mtxResult))
  1702. return FALSE;
  1703. double* x = mtxResult.GetData();
  1704.     
  1705. q=1.0+eps;
  1706.     while (q>=eps)
  1707.     { 
  1708. // 迭代次数已达最大值，仍为求得结果，求解失败
  1709. if (i==0)
  1710. return FALSE;
  1711.         
  1712. // 迭代次数减1
  1713. i=i-1;
  1715. // 矩阵运算
  1716. CMatrix mtxE = m_mtxCoef*mtxResult;
  1717. CMatrix mtxR = m_mtxConst - mtxE;
  1719. // 用全选主元高斯消元法求解
  1720. CLEquations leqs(m_mtxCoef, mtxR);
  1721. CMatrix mtxRR;
  1722. if (! leqs.GetRootsetGauss(mtxRR))
  1723. return FALSE;
  1725. double* r = mtxRR.GetData();
  1726.         
  1727. q=0.0;
  1728.         for ( k=0; k<=n-1; k++)
  1729.         { 
  1730. qq=fabs(r[k])/(1.0+fabs(x[k]+r[k]));
  1731.             if (qq>q) 
  1732. q=qq;
  1733.         }
  1734.         
  1735. for ( k=0; k<=n-1; k++) 
  1736. x[k]=x[k]+r[k];
  1738.     }
  1739.     
  1740. // 求解成功
  1741. return TRUE;
  1742. }