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

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

数学计算

开发平台：

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. CLEquations::CLEquations()
  23. {
  24. }
  26. CLEquations::CLEquations(const CMatrix& mtxCoef, const CMatrix& mtxConst)
  27. {
  28. ASSERT(Init(mtxCoef, mtxConst));
  29. }
  31. CLEquations::~CLEquations()
  32. {
  33. }
  35. BOOL CLEquations::Init(const CMatrix& mtxCoef, const CMatrix& mtxConst)
  36. {
  37. if (mtxCoef.GetNumRows() != mtxConst.GetNumRows())
  38. return FALSE;
  40. m_mtxCoef = mtxCoef;
  41. m_mtxConst = mtxConst;
  43. return TRUE;
  44. }
  46. // 全选主元高斯消去法
  47. BOOL CLEquations::GetRootsetGauss(CMatrix& mtxResult)
  48. { 
  49. int *pnJs,l,k,i,j,nIs,p,q;
  50.     double d,t;
  52. mtxResult = m_mtxConst;
  53. double *pDataCoef = m_mtxCoef.GetData();
  54. double *pDataConst = mtxResult.GetData();
  55. int n = GetNumUnknowns();
  57.     pnJs = new int[n];
  59.     l=1;
  60.     for (k=0;k<=n-2;k++)
  61.     { 
  62. d=0.0;
  63.         for (i=k;i<=n-1;i++)
  64. {
  65. for (j=k;j<=n-1;j++)
  66.             { 
  67. t=fabs(pDataCoef[i*n+j]);
  68. if (t>d) 
  69. { 
  70. d=t; 
  71. pnJs[k]=j; 
  72. nIs=i;
  73. }
  74.             }
  75. }
  77.         if (d == 0.0) 
  78. l=0;
  79.         else
  80.         { 
  81. if (pnJs[k]!=k)
  82.             {
  83. for (i=0;i<=n-1;i++)
  84.                 { 
  85. p=i*n+k; 
  86. q=i*n+pnJs[k];
  87. t=pDataCoef[p]; 
  88. pDataCoef[p]=pDataCoef[q]; 
  89. pDataCoef[q]=t;
  90.                 }
  91. }
  93.             if (nIs!=k)
  94.             { 
  95. for (j=k;j<=n-1;j++)
  96.                 { 
  97. p=k*n+j; 
  98. q=nIs*n+j;
  99.                     t=pDataCoef[p]; 
  100. pDataCoef[p]=pDataCoef[q]; 
  101. pDataCoef[q]=t;
  102.                 }
  103.                 
  104. t=pDataConst[k]; 
  105. pDataConst[k]=pDataConst[nIs]; 
  106. pDataConst[nIs]=t;
  107.             }
  108.         }
  109.         
  110. if (l==0)
  111.         { 
  112. delete[] pnJs;
  113.             return FALSE;
  114.         }
  115.         
  116. d=pDataCoef[k*n+k];
  117.         for (j=k+1;j<=n-1;j++)
  118.         { 
  119. p=k*n+j; 
  120. pDataCoef[p]=pDataCoef[p]/d;
  121. }
  122.         
  123. pDataConst[k]=pDataConst[k]/d;
  124.         for (i=k+1;i<=n-1;i++)
  125.         { 
  126. for (j=k+1;j<=n-1;j++)
  127.             { 
  128. p=i*n+j;
  129.                 pDataCoef[p]=pDataCoef[p]-pDataCoef[i*n+k]*pDataCoef[k*n+j];
  130.             }
  131.             
  132. pDataConst[i]=pDataConst[i]-pDataCoef[i*n+k]*pDataConst[k];
  133.         }
  134.     }
  135.     
  136. d=pDataCoef[(n-1)*n+n-1];
  137.     if (d == 0.0)
  138. { 
  139. delete[] pnJs;
  140. return FALSE;
  141. }
  143.     pDataConst[n-1]=pDataConst[n-1]/d;
  144.     for (i=n-2;i>=0;i--)
  145.     { 
  146. t=0.0;
  147.         for (j=i+1;j<=n-1;j++)
  148. t=t+pDataCoef[i*n+j]*pDataConst[j];
  149.         pDataConst[i]=pDataConst[i]-t;
  150.     }
  151.     
  152. pnJs[n-1]=n-1;
  153.     for (k=n-1;k>=0;k--)
  154.     {
  155. if (pnJs[k]!=k)
  156.         { 
  157. t=pDataConst[k]; 
  158. pDataConst[k]=pDataConst[pnJs[k]]; 
  159. pDataConst[pnJs[k]]=t;
  160. }
  161. }
  163.     delete[] pnJs;
  165.     return TRUE;
  166. }
  168. // 全选主元高斯－约当消去法
  169. BOOL CLEquations::GetRootsetGaussJordan(CMatrix& mtxResult)
  170. { 
  171. int *pnJs,l,k,i,j,nIs,p,q;
  172.     double d,t;
  174. mtxResult = m_mtxConst;
  175. double *pDataCoef = m_mtxCoef.GetData();
  176. double *pDataConst = mtxResult.GetData();
  177. int n = GetNumUnknowns();
  178. int m = m_mtxConst.GetNumColumns();
  180.     pnJs = new int[n];
  182.     l=1;
  183.     for (k=0;k<=n-1;k++)
  184.     { 
  185. d=0.0;
  186.         for (i=k;i<=n-1;i++)
  187.         {
  188. for (j=k;j<=n-1;j++)
  189.             { 
  190. t=fabs(pDataCoef[i*n+j]);
  191. if (t>d) 
  192. { 
  193. d=t; 
  194. pnJs[k]=j; 
  195. nIs=i;
  196. }
  197.             }
  198. }
  200.         if (d+1.0==1.0) 
  201. l=0;
  202.         else
  203.         { 
  204. if (pnJs[k]!=k)
  205.             {
  206. for (i=0;i<=n-1;i++)
  207.                 { 
  208. p=i*n+k; 
  209. q=i*n+pnJs[k];
  210. t=pDataCoef[p]; 
  211. pDataCoef[p]=pDataCoef[q]; 
  212. pDataCoef[q]=t;
  213.                 }
  214. }
  216.             if (nIs!=k)
  217.             { 
  218. for (j=k;j<=n-1;j++)
  219.                 { 
  220. p=k*n+j; 
  221. q=nIs*n+j;
  222.                     t=pDataCoef[p]; 
  223. pDataCoef[p]=pDataCoef[q]; 
  224. pDataCoef[q]=t;
  225.                 }
  226.                 
  227. for (j=0;j<=m-1;j++)
  228.                 { 
  229. p=k*m+j; 
  230. q=nIs*m+j;
  231.                     t=pDataConst[p]; 
  232. pDataConst[p]=pDataConst[q]; 
  233. pDataConst[q]=t;
  234.                 }
  235.             }
  236.         }
  237.         
  238. if (l==0)
  239.         { 
  240. delete[] pnJs;
  241.             return FALSE;
  242.         }
  243.         
  244. d=pDataCoef[k*n+k];
  245.         for (j=k+1;j<=n-1;j++)
  246.         { 
  247. p=k*n+j; 
  248. pDataCoef[p]=pDataCoef[p]/d;
  249. }
  250.         
  251. for (j=0;j<=m-1;j++)
  252.         { 
  253. p=k*m+j; 
  254. pDataConst[p]=pDataConst[p]/d;
  255. }
  256.         
  257. for (j=k+1;j<=n-1;j++)
  258.         {
  259. for (i=0;i<=n-1;i++)
  260.             { 
  261. p=i*n+j;
  262. if (i!=k)
  263. pDataCoef[p]=pDataCoef[p]-pDataCoef[i*n+k]*pDataCoef[k*n+j];
  264.             }
  265. }
  267.         for (j=0;j<=m-1;j++)
  268.         {
  269. for (i=0;i<=n-1;i++)
  270. { 
  271. p=i*m+j;
  272. if (i!=k)
  273. pDataConst[p]=pDataConst[p]-pDataCoef[i*n+k]*pDataConst[k*m+j];
  274. }
  275. }
  276.     }
  277.     
  278. for (k=n-1;k>=0;k--)
  279.     {
  280. if (pnJs[k]!=k)
  281. {
  282. for (j=0;j<=m-1;j++)
  283. { 
  284. p=k*m+j; 
  285. q=pnJs[k]*m+j;
  286. t=pDataConst[p]; 
  287. pDataConst[p]=pDataConst[q]; 
  288. pDataConst[q]=t;
  289. }
  290. }
  291. }
  293.     delete[] pnJs;
  294.     
  295. return TRUE;
  296. }
  298. // 复系数方程组的全选主元高斯消去法
  299. BOOL CLEquations::GetRootsetGauss(CMatrix mtxCoefImag, CMatrix mtxConstImag, CMatrix& mtxResult, CMatrix& mtxResultImag)
  300. { 
  301. int *pnJs,l,k,i,j,nIs,u,v;
  302.     double p,q,s,d;
  304. mtxResult = m_mtxConst;
  305. mtxResultImag = mtxConstImag;
  306. double *pDataCoef = m_mtxCoef.GetData();
  307. double *pDataConst = mtxResult.GetData();
  308. double *pDataCoefImag = mtxCoefImag.GetData();
  309. double *pDataConstImag = mtxResultImag.GetData();
  310. int n = GetNumUnknowns();
  311. int m = m_mtxConst.GetNumColumns();
  313.     pnJs = new int[n];
  314.     
  315. for (k=0;k<=n-2;k++)
  316.     { 
  317. d=0.0;
  318.         for (i=k;i<=n-1;i++)
  319.         {
  320. for (j=k;j<=n-1;j++)
  321. { 
  322. u=i*n+j;
  323. p=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
  324. if (p>d) 
  325. {
  326. d=p;
  327. pnJs[k]=j;
  328. nIs=i;
  329. }
  330. }
  331.         }
  332.         
  333. if (d == 0.0)
  334. {
  335. delete[] pnJs;
  336. return FALSE;
  337.         }
  338.         
  339. if (nIs!=k)
  340.         { 
  341. for (j=k;j<=n-1;j++)
  342.             { 
  343. u=k*n+j; 
  344. v=nIs*n+j;
  345.                 p=pDataCoef[u]; 
  346. pDataCoef[u]=pDataCoef[v]; 
  347. pDataCoef[v]=p;
  348.                 p=pDataCoefImag[u]; 
  349. pDataCoefImag[u]=pDataCoefImag[v]; 
  350. pDataCoefImag[v]=p;
  351.             }
  352.             
  353. p=pDataConst[k]; 
  354. pDataConst[k]=pDataConst[nIs]; 
  355. pDataConst[nIs]=p;
  356.             p=pDataConstImag[k]; 
  357. pDataConstImag[k]=pDataConstImag[nIs]; 
  358. pDataConstImag[nIs]=p;
  359.         }
  360.         
  361. if (pnJs[k]!=k)
  362.         {
  363. for (i=0;i<=n-1;i++)
  364.             { 
  365. u=i*n+k; 
  366. v=i*n+pnJs[k];
  367. p=pDataCoef[u]; 
  368. pDataCoef[u]=pDataCoef[v]; 
  369. pDataCoef[v]=p;
  370. p=pDataCoefImag[u]; 
  371. pDataCoefImag[u]=pDataCoefImag[v]; 
  372. pDataCoefImag[v]=p;
  373.             }
  374. }
  376.         v=k*n+k;
  377.         for (j=k+1;j<=n-1;j++)
  378.         { 
  379. u=k*n+j;
  380.             p=pDataCoef[u]*pDataCoef[v]; 
  381. q=-pDataCoefImag[u]*pDataCoefImag[v];
  382.             s=(pDataCoef[v]-pDataCoefImag[v])*(pDataCoef[u]+pDataCoefImag[u]);
  383.             pDataCoef[u]=(p-q)/d; 
  384. pDataCoefImag[u]=(s-p-q)/d;
  385.         }
  386.         
  387. p=pDataConst[k]*pDataCoef[v]; 
  388. q=-pDataConstImag[k]*pDataCoefImag[v];
  389.         s=(pDataCoef[v]-pDataCoefImag[v])*(pDataConst[k]+pDataConstImag[k]);
  390.         pDataConst[k]=(p-q)/d; 
  391. pDataConstImag[k]=(s-p-q)/d;
  393.         for (i=k+1;i<=n-1;i++)
  394.         { 
  395. u=i*n+k;
  396.             for (j=k+1;j<=n-1;j++)
  397.             { 
  398. v=k*n+j; 
  399. l=i*n+j;
  400.                 p=pDataCoef[u]*pDataCoef[v]; 
  401. q=pDataCoefImag[u]*pDataCoefImag[v];
  402.                 s=(pDataCoef[u]+pDataCoefImag[u])*(pDataCoef[v]+pDataCoefImag[v]);
  403.                 pDataCoef[l]=pDataCoef[l]-p+q;
  404.                 pDataCoefImag[l]=pDataCoefImag[l]-s+p+q;
  405.             }
  406.             
  407. p=pDataCoef[u]*pDataConst[k]; 
  408. q=pDataCoefImag[u]*pDataConstImag[k];
  409.             s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[k]+pDataConstImag[k]);
  410.             pDataConst[i]=pDataConst[i]-p+q; 
  411. pDataConstImag[i]=pDataConstImag[i]-s+p+q;
  412.         }
  413.     }
  414.     
  415. u=(n-1)*n+n-1;
  416.     d=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
  418.     if (d == 0.0)
  419. {
  420. delete[] pnJs;
  421. return FALSE;
  422.     }
  424.     p=pDataCoef[u]*pDataConst[n-1]; q=-pDataCoefImag[u]*pDataConstImag[n-1];
  425.     s=(pDataCoef[u]-pDataCoefImag[u])*(pDataConst[n-1]+pDataConstImag[n-1]);
  426.     pDataConst[n-1]=(p-q)/d; pDataConstImag[n-1]=(s-p-q)/d;
  428.     for (i=n-2;i>=0;i--)
  429.     {
  430. for (j=i+1;j<=n-1;j++)
  431. { 
  432. u=i*n+j;
  433. p=pDataCoef[u]*pDataConst[j]; 
  434. q=pDataCoefImag[u]*pDataConstImag[j];
  435. s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[j]+pDataConstImag[j]);
  436. pDataConst[i]=pDataConst[i]-p+q;
  437. pDataConstImag[i]=pDataConstImag[i]-s+p+q;
  438. }
  439.     }
  441.     pnJs[n-1]=n-1;
  442.     for (k=n-1;k>=0;k--)
  443.     {
  444. if (pnJs[k]!=k)
  445.         { 
  446. p=pDataConst[k]; 
  447. pDataConst[k]=pDataConst[pnJs[k]]; 
  448. pDataConst[pnJs[k]]=p;
  449. p=pDataConstImag[k]; 
  450. pDataConstImag[k]=pDataConstImag[pnJs[k]]; 
  451. pDataConstImag[pnJs[k]]=p;
  452.         }
  453. }
  455. delete[] pnJs;
  457. return TRUE;
  458. }
  460. // 复系数方程组的全选主元高斯－约当消去法
  461. BOOL CLEquations::GetRootsetGaussJordan(CMatrix mtxCoefImag, CMatrix mtxConstImag, CMatrix& mtxResult, CMatrix& mtxResultImag)
  462. {
  463. int *pnJs,l,k,i,j,nIs,u,v;
  464.     double p,q,s,d;
  466. mtxResult = m_mtxConst;
  467. mtxResultImag = mtxConstImag;
  468. double *pDataCoef = m_mtxCoef.GetData();
  469. double *pDataConst = mtxResult.GetData();
  470. double *pDataCoefImag = mtxCoefImag.GetData();
  471. double *pDataConstImag = mtxResultImag.GetData();
  472. int n = GetNumUnknowns();
  473. int m = m_mtxConst.GetNumColumns();
  475.     pnJs = new int[n];
  477.     for (k=0;k<=n-1;k++)
  478.     { 
  479. d=0.0;
  480.         for (i=k;i<=n-1;i++)
  481.         {
  482. for (j=k;j<=n-1;j++)
  483. { 
  484. u=i*n+j;
  485. p=pDataCoef[u]*pDataCoef[u]+pDataCoefImag[u]*pDataCoefImag[u];
  486. if (p>d) 
  487. {
  488. d=p;
  489. pnJs[k]=j;
  490. nIs=i;
  491. }
  492. }
  493.         }
  494.         
  495. if (d == 0.0)
  496.         {
  497. delete[] pnJs;
  498.             return FALSE;
  499.         }
  500.         
  501. if (nIs!=k)
  502.         { 
  503. for (j=k;j<=n-1;j++)
  504.             { 
  505. u=k*n+j; 
  506. v=nIs*n+j;
  507.                 p=pDataCoef[u]; 
  508. pDataCoef[u]=pDataCoef[v]; 
  509. pDataCoef[v]=p;
  510.                 p=pDataCoefImag[u]; 
  511. pDataCoefImag[u]=pDataCoefImag[v]; 
  512. pDataCoefImag[v]=p;
  513.             }
  514.             
  515. for (j=0;j<=m-1;j++)
  516.             { 
  517. u=k*m+j; 
  518. v=nIs*m+j;
  519.                 p=pDataConst[u]; 
  520. pDataConst[u]=pDataConst[v]; 
  521. pDataConst[v]=p;
  522.                 p=pDataConstImag[u]; 
  523. pDataConstImag[u]=pDataConstImag[v]; 
  524. pDataConstImag[v]=p;
  525.             }
  526.         }
  527.         
  528. if (pnJs[k]!=k)
  529.         {
  530. for (i=0;i<=n-1;i++)
  531.             { 
  532. u=i*n+k; 
  533. v=i*n+pnJs[k];
  534. p=pDataCoef[u]; 
  535. pDataCoef[u]=pDataCoef[v]; 
  536. pDataCoef[v]=p;
  537. p=pDataCoefImag[u]; 
  538. pDataCoefImag[u]=pDataCoefImag[v]; 
  539. pDataCoefImag[v]=p;
  540.             }
  541. }
  543.         v=k*n+k;
  544.         for (j=k+1;j<=n-1;j++)
  545.         { 
  546. u=k*n+j;
  547.             p=pDataCoef[u]*pDataCoef[v]; 
  548. q=-pDataCoefImag[u]*pDataCoefImag[v];
  549.             s=(pDataCoef[v]-pDataCoefImag[v])*(pDataCoef[u]+pDataCoefImag[u]);
  550.             pDataCoef[u]=(p-q)/d; 
  551. pDataCoefImag[u]=(s-p-q)/d;
  552.         }
  553.         
  554. for (j=0;j<=m-1;j++)
  555.         { 
  556. u=k*m+j;
  557.             p=pDataConst[u]*pDataCoef[v]; 
  558. q=-pDataConstImag[u]*pDataCoefImag[v];
  559.             s=(pDataCoef[v]-pDataCoefImag[v])*(pDataConst[u]+pDataConstImag[u]);
  560.             pDataConst[u]=(p-q)/d; 
  561. pDataConstImag[u]=(s-p-q)/d;
  562.         }
  563.         
  564. for (i=0;i<=n-1;i++)
  565.         {
  566. if (i!=k)
  567.             { 
  568. u=i*n+k;
  569. for (j=k+1;j<=n-1;j++)
  570.                 { 
  571. v=k*n+j; 
  572. l=i*n+j;
  573. p=pDataCoef[u]*pDataCoef[v]; 
  574. q=pDataCoefImag[u]*pDataCoefImag[v];
  575. s=(pDataCoef[u]+pDataCoefImag[u])*(pDataCoef[v]+pDataCoefImag[v]);
  576. pDataCoef[l]=pDataCoef[l]-p+q;
  577. pDataCoefImag[l]=pDataCoefImag[l]-s+p+q;
  578.                 }
  579.             
  580. for (j=0;j<=m-1;j++)
  581. { 
  582. l=i*m+j; 
  583. v=k*m+j;
  584. p=pDataCoef[u]*pDataConst[v]; q=pDataCoefImag[u]*pDataConstImag[v];
  585. s=(pDataCoef[u]+pDataCoefImag[u])*(pDataConst[v]+pDataConstImag[v]);
  586. pDataConst[l]=pDataConst[l]-p+q; 
  587. pDataConstImag[l]=pDataConstImag[l]-s+p+q;
  588. }
  589.             }
  590. }
  591.     }
  593. for (k=n-1;k>=0;k--)
  594.     {
  595. if (pnJs[k]!=k)
  596. {
  597. for (j=0;j<=m-1;j++)
  598. { 
  599. u=k*m+j;
  600. v=pnJs[k]*m+j;
  601. p=pDataConst[u]; 
  602. pDataConst[u]=pDataConst[v]; 
  603. pDataConst[v]=p;
  604. p=pDataConstImag[u]; 
  605. pDataConstImag[u]=pDataConstImag[v]; 
  606. pDataConstImag[v]=p;
  607. }
  608. }
  609. }
  611. delete[] pnJs;
  613. return TRUE;
  614. }
  616. // 求解三对角线方程组的追赶法
  617. BOOL CLEquations::GetRootsetTriDiagonal(CMatrix& mtxResult)
  618. { 
  619. int k,j;
  620.     double s;
  621.     
  622. mtxResult = m_mtxConst;
  623. double *pDataConst = mtxResult.GetData();
  625. int n = GetNumUnknowns();
  626. ASSERT(m_mtxCoef.GetNumRows() == n);
  628. double* pDiagData = new double[3*n-2];
  630. k = j = 0;
  631. if (k == 0)
  632. {
  633. pDiagData[j++] = m_mtxCoef.GetElement(k, k);
  634. pDiagData[j++] = m_mtxCoef.GetElement(k, k+1);
  635. }
  636.     for (k=1; k<n-1; ++k)
  637. {
  638. pDiagData[j++] = m_mtxCoef.GetElement(k, k-1);
  639. pDiagData[j++] = m_mtxCoef.GetElement(k, k);
  640. pDiagData[j++] = m_mtxCoef.GetElement(k, k+1);
  641. }
  642. if (k == n-1)
  643. {
  644. pDiagData[j++] = m_mtxCoef.GetElement(k, k-1);
  645. pDiagData[j++] = m_mtxCoef.GetElement(k, k);
  646. }
  648. for (k=0;k<=n-2;k++)
  649.     { 
  650. j=3*k; 
  651. s=pDiagData[j];
  653.         if (fabs(s)+1.0==1.0)
  654. {
  655. delete[] pDiagData;
  656. return FALSE;
  657. }
  659.         pDiagData[j+1]=pDiagData[j+1]/s;
  660.         pDataConst[k]=pDataConst[k]/s;
  661.         pDiagData[j+3]=pDiagData[j+3]-pDiagData[j+2]*pDiagData[j+1];
  662.         pDataConst[k+1]=pDataConst[k+1]-pDiagData[j+2]*pDataConst[k];
  663.     }
  664.     
  665. s=pDiagData[3*n-3];
  666.     if (s == 0.0)
  667. {
  668. delete[] pDiagData;
  669. return FALSE;
  670. }
  671.     
  672. pDataConst[n-1]=pDataConst[n-1]/s;
  673.     for (k=n-2;k>=0;k--)
  674. pDataConst[k]=pDataConst[k]-pDiagData[3*k+1]*pDataConst[k+1];
  675.     
  676. delete[] pDiagData;
  678. return TRUE;
  679. }
  681. // 一般带型方程组的求解
  682. BOOL CLEquations::GetRootsetBand(int nBandWidth, CMatrix& mtxResult)
  683. { 
  684. int ls,k,i,j,is,u,v;
  685.     double p,t;
  686.     
  687. if ((nBandWidth-1)%2 != 0)
  688. return FALSE;
  690. mtxResult = m_mtxConst;
  691. double *pDataConst = mtxResult.GetData();
  693. int m = m_mtxConst.GetNumColumns();
  694. int n = GetNumUnknowns();
  695. ASSERT(m_mtxCoef.GetNumRows() == n);
  697. double* pBandData = new double[nBandWidth*n];
  698. // 半带宽
  699.     ls = (nBandWidth-1)/2;
  701. // 构造带宽数组
  702.     for (i=0; i<n; ++i)
  703. {
  704. j = 0;
  705. for (k=max(0, i-ls); k<max(0, i-ls)+nBandWidth; ++k)
  706. {
  707. if (k < n)
  708. pBandData[i*nBandWidth+j++] = m_mtxCoef.GetElement(i, k);
  709. else
  710. pBandData[i*nBandWidth+j++] = 0;
  711. }
  712. }
  714.     for (k=0;k<=n-2;k++)
  715.     { 
  716. p=0.0;
  717.         for (i=k;i<=ls;i++)
  718.         { 
  719. t=fabs(pBandData[i*nBandWidth]);
  720.             if (t>p) 
  721. {
  722. p=t; 
  723. is=i;
  724. }
  725.         }
  726.         
  727. if (p == 0.0)
  728. {
  729. delete[] pBandData;
  730. return FALSE;
  731. }
  733.         for (j=0;j<=m-1;j++)
  734.         { 
  735. u=k*m+j; 
  736. v=is*m+j;
  737.             t=pDataConst[u]; 
  738. pDataConst[u]=pDataConst[v]; 
  739. pDataConst[v]=t;
  740.         }
  741.         
  742. for (j=0;j<=nBandWidth-1;j++)
  743.         { 
  744. u=k*nBandWidth+j; 
  745. v=is*nBandWidth+j;
  746.             t=pBandData[u]; 
  747. pBandData[u]=pBandData[v]; 
  748. pBandData[v]=t;
  749.         }
  750.         
  751. for (j=0;j<=m-1;j++)
  752.         { 
  753. u=k*m+j; 
  754. pDataConst[u]=pDataConst[u]/pBandData[k*nBandWidth];
  755. }
  756.         
  757. for (j=1;j<=nBandWidth-1;j++)
  758.         { 
  759. u=k*nBandWidth+j; 
  760. pBandData[u]=pBandData[u]/pBandData[k*nBandWidth];
  761. }
  762.         
  763. for (i=k+1;i<=ls;i++)
  764.         { 
  765. t=pBandData[i*nBandWidth];
  766.             for (j=0;j<=m-1;j++)
  767.             { 
  768. u=i*m+j; 
  769. v=k*m+j;
  770.                 pDataConst[u]=pDataConst[u]-t*pDataConst[v];
  771.             }
  772.             
  773. for (j=1;j<=nBandWidth-1;j++)
  774.             { 
  775. u=i*nBandWidth+j; 
  776. v=k*nBandWidth+j;
  777.                 pBandData[u-1]=pBandData[u]-t*pBandData[v];
  778.             }
  779.             
  780. u=i*nBandWidth+nBandWidth-1; pBandData[u]=0.0;
  781.         }
  782.         
  783. if (ls!=(n-1)) 
  784. ls=ls+1;
  785.     }
  786.     
  787. p=pBandData[(n-1)*nBandWidth];
  788.     if (p == 0.0)
  789. {
  790. delete[] pBandData;
  791. return FALSE;
  792. }
  794.     for (j=0;j<=m-1;j++)
  795.     { 
  796. u=(n-1)*m+j; 
  797. pDataConst[u]=pDataConst[u]/p;
  798. }
  799.     
  800. ls=1;
  801.     for (i=n-2;i>=0;i--)
  802.     { 
  803. for (k=0;k<=m-1;k++)
  804.         { 
  805. u=i*m+k;
  806.             for (j=1;j<=ls;j++)
  807.             { 
  808. v=i*nBandWidth+j; 
  809. is=(i+j)*m+k;
  810.                 pDataConst[u]=pDataConst[u]-pBandData[v]*pDataConst[is];
  811.             }
  812.         }
  813.         
  814. if (ls!=(nBandWidth-1)) 
  815. ls=ls+1;
  816.     }
  817.     
  818. delete[] pBandData;
  820. return TRUE;
  821. }
  823. // 求解对称方程组的分解法
  824. BOOL CLEquations::GetRootsetDjn(CMatrix& mtxResult)
  825. { 
  826. int i,j,l,k,u,v,w,k1,k2,k3;
  827.     double p;
  829. CMatrix mtxCoef = m_mtxCoef;
  830. mtxResult = m_mtxConst;
  831. int n = mtxCoef.GetNumColumns();
  832. int m = mtxResult.GetNumColumns();
  833. double* pDataCoef = mtxCoef.GetData();
  834. double* pDataConst = mtxResult.GetData();
  836.     if (pDataCoef[0] == 0.0)
  837. return FALSE;
  839.     for (i=1; i<=n-1; i++)
  840.     { 
  841. u=i*n; 
  842. pDataCoef[u]=pDataCoef[u]/pDataCoef[0];
  843. }
  844.     
  845. for (i=1; i<=n-2; i++)
  846.     { 
  847. u=i*n+i;
  848.         for (j=1; j<=i; j++)
  849.         { 
  850. v=i*n+j-1; 
  851. l=(j-1)*n+j-1;
  852.             pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v]*pDataCoef[l];
  853.         }
  854.         
  855. p=pDataCoef[u];
  856.         if (p == 0.0)
  857. return FALSE;
  859.         for (k=i+1; k<=n-1; k++)
  860.         { 
  861. u=k*n+i;
  862.             for (j=1; j<=i; j++)
  863.             { 
  864. v=k*n+j-1; 
  865. l=i*n+j-1; 
  866. w=(j-1)*n+j-1;
  867. pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[l]*pDataCoef[w];
  868.             }
  869.             
  870. pDataCoef[u]=pDataCoef[u]/p;
  871.         }
  872.     }
  873.     
  874. u=n*n-1;
  875.     for (j=1; j<=n-1; j++)
  876.     { 
  877. v=(n-1)*n+j-1; 
  878. w=(j-1)*n+j-1;
  879.         pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v]*pDataCoef[w];
  880.     }
  881.     
  882. p=pDataCoef[u];
  883.     if (p == 0.0)
  884. return FALSE;
  886.     for (j=0; j<=m-1; j++)
  887. {
  888.     for (i=1; i<=n-1; i++)
  889. { 
  890. u=i*m+j;
  891. for (k=1; k<=i; k++)
  892. { 
  893. v=i*n+k-1; 
  894. w=(k-1)*m+j;
  895. pDataConst[u]=pDataConst[u]-pDataCoef[v]*pDataConst[w];
  896. }
  897. }
  898. }
  900.     for (i=1; i<=n-1; i++)
  901.     { 
  902. u=(i-1)*n+i-1;
  903.         for (j=i; j<=n-1; j++)
  904.         { 
  905. v=(i-1)*n+j; 
  906. w=j*n+i-1;
  907.             pDataCoef[v]=pDataCoef[u]*pDataCoef[w];
  908.         }
  909.     }
  910.     
  911. for (j=0; j<=m-1; j++)
  912.     { 
  913. u=(n-1)*m+j;
  914.         pDataConst[u]=pDataConst[u]/p;
  915.         for (k=1; k<=n-1; k++)
  916.         { 
  917. k1=n-k; 
  918. k3=k1-1; 
  919. u=k3*m+j;
  920.             for (k2=k1; k2<=n-1; k2++)
  921.             { 
  922. v=k3*n+k2; 
  923. w=k2*m+j;
  924.                 pDataConst[u]=pDataConst[u]-pDataCoef[v]*pDataConst[w];
  925.             }
  926.             
  927. pDataConst[u]=pDataConst[u]/pDataCoef[k3*n+k3];
  928.         }
  929.     }
  930.     
  931. return TRUE;
  932. }
  934. // 求解对称正定方程组的平方根法
  935. BOOL CLEquations::GetRootsetCholesky(CMatrix& mtxResult)
  936. { 
  937. int i,j,k,u,v;
  938.     
  939. CMatrix mtxCoef = m_mtxCoef;
  940. mtxResult = m_mtxConst;
  941. int n = mtxCoef.GetNumColumns();
  942. int m = mtxResult.GetNumColumns();
  943. double* pDataCoef = mtxCoef.GetData();
  944. double* pDataConst = mtxResult.GetData();
  945.     
  946. if (pDataCoef[0] <= 0.0)
  947. return FALSE;
  949. pDataCoef[0]=sqrt(pDataCoef[0]);
  950.     for (j=1; j<=n-1; j++) 
  951. pDataCoef[j]=pDataCoef[j]/pDataCoef[0];
  952.     
  953. for (i=1; i<=n-1; i++)
  954.     { 
  955. u=i*n+i;
  956.         for (j=1; j<=i; j++)
  957.         { 
  958. v=(j-1)*n+i;
  959.             pDataCoef[u]=pDataCoef[u]-pDataCoef[v]*pDataCoef[v];
  960.         }
  961.         
  962. if (pDataCoef[u] <= 0.0)
  963. return FALSE;
  965.         pDataCoef[u]=sqrt(pDataCoef[u]);
  966.         if (i!=(n-1))
  967.         { 
  968. for (j=i+1; j<=n-1; j++)
  969.             { 
  970. v=i*n+j;
  971.                 for (k=1; k<=i; k++)
  972. pDataCoef[v]=pDataCoef[v]-pDataCoef[(k-1)*n+i]*pDataCoef[(k-1)*n+j];
  973.                 pDataCoef[v]=pDataCoef[v]/pDataCoef[u];
  974.             }
  975.         }
  976.     }
  977.     
  978. for (j=0; j<=m-1; j++)
  979.     { 
  980. pDataConst[j]=pDataConst[j]/pDataCoef[0];
  981.         for (i=1; i<=n-1; i++)
  982.         { 
  983. u=i*n+i; 
  984. v=i*m+j;
  985.             for (k=1; k<=i; k++)
  986. pDataConst[v]=pDataConst[v]-pDataCoef[(k-1)*n+i]*pDataConst[(k-1)*m+j];
  987.             pDataConst[v]=pDataConst[v]/pDataCoef[u];
  988.         }
  989.     }
  990.     
  991. for (j=0; j<=m-1; j++)
  992.     { 
  993. u=(n-1)*m+j;
  994.         pDataConst[u]=pDataConst[u]/pDataCoef[n*n-1];
  995.         for (k=n-1; k>=1; k--)
  996.         { 
  997. u=(k-1)*m+j;
  998.             for (i=k; i<=n-1; i++)
  999.             { 
  1000. v=(k-1)*n+i;
  1001.                 pDataConst[u]=pDataConst[u]-pDataCoef[v]*pDataConst[i*m+j];
  1002.             }
  1003.             
  1004. v=(k-1)*n+k-1;
  1005.             pDataConst[u]=pDataConst[u]/pDataCoef[v];
  1006.         }
  1007.     }
  1008.     
  1009. return TRUE;
  1010. }
  1012. // 求解大型稀疏方程组的全选主元高斯－约去消去法
  1013. BOOL CLEquations::GetRootsetGgje(CMatrix& mtxResult)
  1014. { 
  1015. int *pnJs,i,j,k,nIs,u,v;
  1016.     double d,t;
  1017.     
  1018. CMatrix mtxCoef = m_mtxCoef;
  1019. mtxResult = m_mtxConst;
  1020. int n = mtxCoef.GetNumColumns();
  1021. double* pDataCoef = mtxCoef.GetData();
  1022. double* pDataConst = mtxResult.GetData();
  1024. pnJs = new int[n];
  1026.     for (k=0; k<=n-1; k++)
  1027.     { 
  1028. d=0.0;
  1029.         for (i=k; i<=n-1; i++)
  1030. {
  1031. for (j=k; j<=n-1; j++)
  1032. { 
  1033. t=fabs(pDataCoef[i*n+j]);
  1034. if (t>d) 
  1035. {
  1036. d=t; 
  1037. pnJs[k]=j; 
  1038. nIs=i;
  1039. }
  1040. }
  1041. }
  1043.         if (d == 0.0)
  1044. {
  1045. delete[] pnJs;
  1046. return FALSE;
  1047. }
  1049.         if (nIs!=k)
  1050.         { 
  1051. for (j=k; j<=n-1; j++)
  1052.             { 
  1053. u=k*n+j; 
  1054. v=nIs*n+j;
  1055.                 t=pDataCoef[u]; 
  1056. pDataCoef[u]=pDataCoef[v]; 
  1057. pDataCoef[v]=t;
  1058.             }
  1059.             
  1060. t=pDataConst[k]; 
  1061. pDataConst[k]=pDataConst[nIs]; 
  1062. pDataConst[nIs]=t;
  1063.         }
  1064.         
  1065. if (pnJs[k]!=k)
  1066.         {
  1067. for (i=0; i<=n-1; i++)
  1068.             { 
  1069. u=i*n+k; 
  1070. v=i*n+pnJs[k];
  1071. t=pDataCoef[u]; 
  1072. pDataCoef[u]=pDataCoef[v]; 
  1073. pDataCoef[v]=t;
  1074.             }
  1075. }
  1077.         t=pDataCoef[k*n+k];
  1078.         for (j=k+1; j<=n-1; j++)
  1079.         { 
  1080. u=k*n+j;
  1081.             if (pDataCoef[u]!=0.0) 
  1082. pDataCoef[u]=pDataCoef[u]/t;
  1083.         }
  1084.         
  1085. pDataConst[k]=pDataConst[k]/t;
  1086.         for (j=k+1; j<=n-1; j++)
  1087.         { 
  1088. u=k*n+j;
  1089.             if (pDataCoef[u]!=0.0)
  1090.             { 
  1091. for (i=0; i<=n-1; i++)
  1092.                 { 
  1093. v=i*n+k;
  1094.                     if ((i!=k)&&(pDataCoef[v]!=0.0))
  1095.                     { 
  1096. nIs=i*n+j;
  1097.                         pDataCoef[nIs]=pDataCoef[nIs]-pDataCoef[v]*pDataCoef[u];
  1098.                     }
  1099.                 }
  1100.             }
  1101.         }
  1102.         
  1103. for (i=0; i<=n-1; i++)
  1104.         { 
  1105. u=i*n+k;
  1106.             if ((i!=k)&&(pDataCoef[u]!=0.0))
  1107. pDataConst[i]=pDataConst[i]-pDataCoef[u]*pDataConst[k];
  1108.         }
  1109.     }
  1110.     
  1111. for (k=n-1; k>=0; k--)
  1112.     {
  1113. if (k!=pnJs[k])
  1114.         { 
  1115. t=pDataConst[k]; 
  1116. pDataConst[k]=pDataConst[pnJs[k]]; 
  1117. pDataConst[pnJs[k]]=t;
  1118. }
  1119. }
  1121.     delete[] pnJs;
  1122.     
  1123. return TRUE;
  1124. }
  1126. // 求解托伯利兹方程组的列文逊方法
  1127. BOOL CLEquations::GetRootsetTlvs(CMatrix& mtxResult)
  1128. { 
  1129. int i,j,k;
  1130.     double a,beta,q,c,h,*y,*s;
  1132.     int n = m_mtxCoef.GetNumColumns();
  1134. mtxResult.Init(n, 1);
  1135. double* x = mtxResult.GetData();
  1137. double* pDataConst = m_mtxConst.GetData();
  1139. double* t = new double[n];
  1141. // 构造T数组
  1142.     for (i=0; i<n; ++i)
  1143. t[i] = m_mtxCoef.GetElement(0, i);
  1145. s = new double[n];
  1146.     y = new double[n];
  1148.     a=t[0];
  1149.     if (a == 0.0)
  1150.     { 
  1151. delete[] s;
  1152. delete[] y;
  1153. delete[] t;
  1154. return FALSE;
  1155. }
  1157.     y[0]=1.0; 
  1158. x[0]=pDataConst[0]/a;
  1159.     for (k=1; k<=n-1; k++)
  1160.     { 
  1161. beta=0.0; 
  1162. q=0.0;
  1163.         for (j=0; j<=k-1; j++)
  1164.         { 
  1165. beta=beta+y[j]*t[j+1];
  1166.             q=q+x[j]*t[k-j];
  1167.         }
  1168.         
  1169. if (a == 0.0)
  1170. { 
  1171. delete[] s;
  1172. delete[] y;
  1173. delete[] t;
  1174. return FALSE;
  1175. }
  1177.         c=-beta/a; 
  1178. s[0]=c*y[k-1]; 
  1179. y[k]=y[k-1];
  1180.         if (k!=1)
  1181.         {
  1182. for (i=1; i<=k-1; i++)
  1183. s[i]=y[i-1]+c*y[k-i-1];
  1184. }
  1186.         a=a+c*beta;
  1187. if (a == 0.0)
  1188. { 
  1189. delete[] s;
  1190. delete[] y;
  1191. delete[] t;
  1192. return FALSE;
  1193. }
  1195.         h=(pDataConst[k]-q)/a;
  1196.         for (i=0; i<=k-1; i++)
  1197.         { 
  1198. x[i]=x[i]+h*s[i]; 
  1199. y[i]=s[i];
  1200. }
  1201.         
  1202. x[k]=h*y[k];
  1203.     }
  1204.     
  1205. delete[] s;
  1206. delete[] y;
  1207. delete[] t;
  1209.     return TRUE;
  1210. }
  1212. // 高斯－赛德尔迭代法
  1213. BOOL CLEquations::GetRootsetGaussSeidel(CMatrix& mtxResult, double eps /*= 0.000001*/)
  1214. { 
  1215. int i,j,u,v;
  1216.     double p,t,s,q;
  1218.     int n = m_mtxCoef.GetNumColumns();
  1220. mtxResult.Init(n, 1);
  1221. double* x = mtxResult.GetData();
  1223. double* pDataCoef = m_mtxCoef.GetData();
  1224. double* pDataConst = m_mtxConst.GetData();
  1225.     
  1226. for (i=0; i<=n-1; i++)
  1227.     { 
  1228. u=i*n+i; 
  1229. p=0.0; 
  1230. x[i]=0.0;
  1231.         for (j=0; j<=n-1; j++)
  1232.         {
  1233. if (i!=j)
  1234.             { 
  1235. v=i*n+j; 
  1236. p=p+fabs(pDataCoef[v]);
  1237. }
  1238. }
  1240.         if (p>=fabs(pDataCoef[u]))
  1241. return FALSE;
  1242.     }
  1244.     p=eps+1.0;
  1245.     while (p>=eps)
  1246.     { 
  1247. p=0.0;
  1248.         for (i=0; i<=n-1; i++)
  1249.         { 
  1250. t=x[i]; 
  1251. s=0.0;
  1252.             for (j=0; j<=n-1; j++)
  1253. if (j!=i) 
  1254. s=s+pDataCoef[i*n+j]*x[j];
  1256.             x[i]=(pDataConst[i]-s)/pDataCoef[i*n+i];
  1257.             q=fabs(x[i]-t)/(1.0+fabs(x[i]));
  1258.             if (q>p) 
  1259. p=q;
  1260.         }
  1261.     }
  1262.     
  1263. return TRUE;
  1264. }
  1266. // 求解对称正定方程组的共轭梯度法
  1267. void CLEquations::GetRootsetGrad(CMatrix& mtxResult, double eps /*= 0.000001*/)
  1268. { 
  1269. int i,k;
  1270.     double *p,*r,*s,*q,alpha,beta,d,e;
  1272.     int n = GetNumUnknowns();
  1274. mtxResult.Init(n, 1);
  1275. double* x = mtxResult.GetData();
  1277. CMatrix mtxP(n, 1);
  1278.     p = mtxP.GetData();
  1280. double* pDataCoef = m_mtxCoef.GetData();
  1281. double* pDataConst = m_mtxConst.GetData();
  1283.     r = new double[n];
  1285.     for (i=0; i<=n-1; i++)
  1286.     { 
  1287. x[i]=0.0; 
  1288. p[i]=pDataConst[i]; 
  1289. r[i]=pDataConst[i]; 
  1290. }
  1291.     
  1292. i=0;
  1293.     while (i<=n-1)
  1294.     { 
  1295. CMatrix mtxS = m_mtxCoef * mtxP;
  1296.     s = mtxS.GetData();
  1297.         
  1298. d=0.0; 
  1299. e=0.0;
  1300.         for (k=0; k<=n-1; k++)
  1301.         { 
  1302. d=d+p[k]*pDataConst[k]; 
  1303. e=e+p[k]*s[k]; 
  1304. }
  1305.         
  1306. alpha=d/e;
  1307.         for (k=0; k<=n-1; k++)
  1308. x[k]=x[k]+alpha*p[k];
  1309.         
  1310. CMatrix mtxQ = m_mtxCoef * mtxResult;
  1311.     q = mtxQ.GetData();
  1312.         
  1313. d=0.0;
  1314.         for (k=0; k<=n-1; k++)
  1315.         { 
  1316. r[k]=pDataConst[k]-q[k]; 
  1317. d=d+r[k]*s[k]; 
  1318. }
  1319.         
  1320. beta=d/e; d=0.0;
  1321.         for (k=0; k<=n-1; k++) 
  1322. d=d+r[k]*r[k];
  1323.         
  1324. d=sqrt(d);
  1325.         if (d<eps)
  1326. break;
  1328.         for (k=0; k<=n-1; k++)
  1329. p[k]=r[k]-beta*p[k];
  1330.         
  1331. i=i+1;
  1332.     }
  1334. delete[] r;
  1335. }
  1337. // 求解线性最小二乘问题的豪斯荷尔德变换法
  1338. BOOL CLEquations::GetRootsetMqr(CMatrix& mtxResult, CMatrix& mtxQ, CMatrix& mtxR)
  1339. { 
  1340. int i,j;
  1341.     double d;
  1343.     int m = m_mtxCoef.GetNumRows();
  1344.     int n = m_mtxCoef.GetNumColumns();
  1345. if (m < n)
  1346. return FALSE;
  1348. mtxResult = m_mtxConst;
  1349. double* pDataConst = mtxResult.GetData();
  1351. mtxR = m_mtxCoef;
  1352. double* pDataCoef = mtxR.GetData();
  1354.     if (! mtxR.SplitQR(mtxQ))
  1355. return FALSE;
  1357. double* c = new double[n];
  1358. double* q = mtxQ.GetData();
  1360.     for (i=0; i<=n-1; i++)
  1361.     { 
  1362. d=0.0;
  1363.         for (j=0; j<=m-1; j++)
  1364. d=d+q[j*m+i]*pDataConst[j];
  1365.     
  1366. c[i]=d;
  1367.     }
  1368.     
  1369. pDataConst[n-1]=c[n-1]/pDataCoef[n*n-1];
  1370.     for (i=n-2; i>=0; i--)
  1371.     { 
  1372. d=0.0;
  1373.         for (j=i+1; j<=n-1; j++)
  1374. d=d+pDataCoef[i*n+j]*pDataConst[j];
  1375.         
  1376. pDataConst[i]=(c[i]-d)/pDataCoef[i*n+i];
  1377.     }
  1378.     
  1379. delete[] c;
  1381. return TRUE;
  1382. }
  1384. // 求解线性最小二乘问题的广义逆法
  1385. BOOL CLEquations::GetRootsetGinv(CMatrix& mtxResult, CMatrix& mtxAP, CMatrix& mtxU, CMatrix& mtxV, double eps /*= 0.000001*/)
  1386. { 
  1387. int i,j;
  1388.     
  1389.     int m = m_mtxCoef.GetNumRows();
  1390.     int n = m_mtxCoef.GetNumColumns();
  1392. mtxResult.Init(n, 1);
  1394. double* pDataConst = m_mtxConst.GetData();
  1395. double* x = mtxResult.GetData();
  1397. CMatrix mtxA = m_mtxCoef;
  1399.     if (! mtxA.GInvertUV(mtxAP, mtxU, mtxV, eps))
  1400. return FALSE;
  1402. double* pAPData = mtxAP.GetData();
  1404.     for (i=0; i<=n-1; i++)
  1405.     { 
  1406. x[i]=0.0;
  1407.         for (j=0; j<=m-1; j++)
  1408. x[i]=x[i]+pAPData[i*m+j]*pDataConst[j];
  1409.     }
  1410.     
  1411. return TRUE;
  1412. }
  1414. // 病态方程组的求解
  1415. BOOL CLEquations::GetRootsetMorbid(CMatrix& mtxResult, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  1416. { 
  1417. int i, k;
  1418.     double q, qq;
  1419.     
  1420. // 方程的阶数
  1421.     int n = GetNumUnknowns();
  1423. // 设定迭代次数, 缺省为60
  1424. i = nMaxIt;
  1425.     
  1426. // 用全选主元高斯消元法求解
  1427. CLEquations leqs(m_mtxCoef, m_mtxConst);
  1428. if (! leqs.GetRootsetGauss(mtxResult))
  1429. return FALSE;
  1430. double* x = mtxResult.GetData();
  1431.     
  1432. q=1.0+eps;
  1433.     while (q>=eps)
  1434.     { 
  1435. // 迭代次数已达最大值，仍为求得结果，求解失败
  1436. if (i==0)
  1437. return FALSE;
  1438.         
  1439. // 迭代次数减1
  1440. i=i-1;
  1442. // 矩阵运算
  1443. CMatrix mtxE = m_mtxCoef*mtxResult;
  1444. CMatrix mtxR = m_mtxConst - mtxE;
  1446. // 用全选主元高斯消元法求解
  1447. CLEquations leqs(m_mtxCoef, mtxR);
  1448. CMatrix mtxRR;
  1449. if (! leqs.GetRootsetGauss(mtxRR))
  1450. return FALSE;
  1452. double* r = mtxRR.GetData();
  1453.         
  1454. q=0.0;
  1455.         for ( k=0; k<=n-1; k++)
  1456.         { 
  1457. qq=fabs(r[k])/(1.0+fabs(x[k]+r[k]));
  1458.             if (qq>q) 
  1459. q=qq;
  1460.         }
  1461.         
  1462. for ( k=0; k<=n-1; k++) 
  1463. x[k]=x[k]+r[k];
  1465.     }
  1466.     
  1467. // 求解成功
  1468. return TRUE;
  1469. }