?? multiplydlg.cpp
字號:
// multiplyDlg.cpp : implementation file
//
#include <math.h>
#include "stdafx.h"
#include "multiply.h"
#include "multiplyDlg.h"
#include "TSMatrix.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
/////////////////////////////////////////////////////////////////////////////
// CAboutDlg dialog used for App About
class CAboutDlg : public CDialog
{
public:
CAboutDlg();
// Dialog Data
//{{AFX_DATA(CAboutDlg)
enum { IDD = IDD_ABOUTBOX };
//}}AFX_DATA
// ClassWizard generated virtual function overrides
//{{AFX_VIRTUAL(CAboutDlg)
protected:
virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV support
//}}AFX_VIRTUAL
// Implementation
protected:
//{{AFX_MSG(CAboutDlg)
//}}AFX_MSG
DECLARE_MESSAGE_MAP()
};
CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD)
{
//{{AFX_DATA_INIT(CAboutDlg)
//}}AFX_DATA_INIT
}
void CAboutDlg::DoDataExchange(CDataExchange* pDX)
{
CDialog::DoDataExchange(pDX);
//{{AFX_DATA_MAP(CAboutDlg)
//}}AFX_DATA_MAP
}
BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)
//{{AFX_MSG_MAP(CAboutDlg)
// No message handlers
//}}AFX_MSG_MAP
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CMultiplyDlg dialog
CMultiplyDlg::CMultiplyDlg(CWnd* pParent /*=NULL*/)
: CDialog(CMultiplyDlg::IDD, pParent)
{
//{{AFX_DATA_INIT(CMultiplyDlg)
// NOTE: the ClassWizard will add member initialization here
//}}AFX_DATA_INIT
// Note that LoadIcon does not require a subsequent DestroyIcon in Win32
m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}
void CMultiplyDlg::DoDataExchange(CDataExchange* pDX)
{
CDialog::DoDataExchange(pDX);
//{{AFX_DATA_MAP(CMultiplyDlg)
// NOTE: the ClassWizard will add DDX and DDV calls here
//}}AFX_DATA_MAP
}
BEGIN_MESSAGE_MAP(CMultiplyDlg, CDialog)
//{{AFX_MSG_MAP(CMultiplyDlg)
ON_WM_SYSCOMMAND()
ON_WM_PAINT()
ON_WM_QUERYDRAGICON()
//}}AFX_MSG_MAP
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CMultiplyDlg message handlers
BOOL CMultiplyDlg::OnInitDialog()
{
CDialog::OnInitDialog();
// Add "About..." menu item to system menu.
// IDM_ABOUTBOX must be in the system command range.
ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);
ASSERT(IDM_ABOUTBOX < 0xF000);
CMenu* pSysMenu = GetSystemMenu(FALSE);
if (pSysMenu != NULL)
{
CString strAboutMenu;
strAboutMenu.LoadString(IDS_ABOUTBOX);
if (!strAboutMenu.IsEmpty())
{
pSysMenu->AppendMenu(MF_SEPARATOR);
pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
}
}
// Set the icon for this dialog. The framework does this automatically
// when the application's main window is not a dialog
SetIcon(m_hIcon, TRUE); // Set big icon
SetIcon(m_hIcon, FALSE); // Set small icon
// TODO: Add extra initialization here
return TRUE; // return TRUE unless you set the focus to a control
}
void CMultiplyDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
if ((nID & 0xFFF0) == IDM_ABOUTBOX)
{
CAboutDlg dlgAbout;
dlgAbout.DoModal();
}
else
{
CDialog::OnSysCommand(nID, lParam);
}
}
// If you add a minimize button to your dialog, you will need the code below
// to draw the icon. For MFC applications using the document/view model,
// this is automatically done for you by the framework.
void CMultiplyDlg::OnPaint()
{
if (IsIconic())
{
CPaintDC dc(this); // device context for painting
SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);
// Center icon in client rectangle
int cxIcon = GetSystemMetrics(SM_CXICON);
int cyIcon = GetSystemMetrics(SM_CYICON);
CRect rect;
GetClientRect(&rect);
int x = (rect.Width() - cxIcon + 1) / 2;
int y = (rect.Height() - cyIcon + 1) / 2;
// Draw the icon
dc.DrawIcon(x, y, m_hIcon);
}
else
{
CDialog::OnPaint();
}
}
// The system calls this to obtain the cursor to display while the user drags
// the minimized window.
HCURSOR CMultiplyDlg::OnQueryDragIcon()
{
return (HCURSOR) m_hIcon;
}
void CMultiplyDlg::OnOK()
{
// TODO: Add extra validation here
U0 = new float[4];
b= new float[4];
ICL.mu = ICL.nu = 3;
ICL.tu = 6;
ICL.data = new Triple[ICL.tu+1];
// ICL.data[1].e = 1.7321;
// ICL.data[1].ii = 1;
// ICL.data[1].jj = 1;
//
// ICL.data[2].e = 0.5774;
// ICL.data[2].ii = 2;
// ICL.data[2].jj = 1;
//
// ICL.data[3].e = 1.2910;
// ICL.data[3].ii = 2;
// ICL.data[3].jj = 2;
//
// ICL.data[4].e = 1.1547;
// ICL.data[4].ii = 3;
// ICL.data[4].jj = 1;
//
// ICL.data[5].e = 0.7746;
// ICL.data[5].ii = 3;
// ICL.data[5].jj = 2;
//
// ICL.data[6].e = 1.0328;
// ICL.data[6].ii = 3;
// ICL.data[6].jj = 3;
U0[1] = 0;
U0[2] = 0;
U0[3] = 0;
b[1] = 3;
b[2] = 2;
b[3] = 5;
K.mu = K.nu = 3;
K.tu = 9;
K.data = new Triple[K.tu];
K.data[1].e = 3;
K.data[1].ii = 1;
K.data[1].jj = 1;
K.data[2].e =1;
K.data[2].ii = 1;
K.data[2].jj = 2;
K.data[3].e =2;
K.data[3].ii = 1;
K.data[3].jj = 3;
K.data[4].e = 1;
K.data[4].ii = 2;
K.data[4].jj = 1;
K.data[5].e = 2;
K.data[5].ii = 2;
K.data[5].jj = 2;
K.data[6].e =1 ;
K.data[6].ii = 2;
K.data[6].jj = 3;
K.data[7].e =2;
K.data[7].ii = 3;
K.data[7].jj = 1;
K.data[8].e =1;
K.data[8].ii = 3;
K.data[8].jj = 2;
K.data[9].e =3;
K.data[9].ii = 3;
K.data[9].jj = 3;
iCholesky();
ICPCG();
CDialog::OnOK();
}
void CMultiplyDlg::offdiagsolve(TSMatrix &M, float *r, float *solv)
{
int mrow;
int i,tp,t;
int* mnum = new int[M.mu+1];
int* mpot = new int[M.mu+1];
for(mrow=1;mrow<=M.mu;++mrow)
mnum[mrow] = 0;
for( t=1;t<=M.tu;++t)
++mnum[M.data[t].ii];
mpot[1] = 1;
//*(cpot+1)=1; //M中第一行第一個非零元在M中的序號為1
for(mrow=2;mrow<=M.mu;++mrow)
mpot[mrow] = mpot[mrow-1] + mnum[mrow-1];
//float si;
float subr;
for (i=1;i<=M.mu;i++)
{
subr = r[i];
if(i<M.mu)
tp=mpot[i+1];
else
tp=M.tu+1;
for(int p=mpot[i];p<tp-1;++p)
{
t = M.data[p].jj;
subr = subr - M.data[p].e * solv[t];
}
if (M.data[tp-1].e==0)
{
solv[i] = 0;
}
else
solv[i] = subr/M.data[tp-1].e;
}
}
void CMultiplyDlg::uppdiagsolve(TSMatrix &M, float *r, float *solv)
{
int mrow;
int i,tp,t;
int* mnum = new int[M.mu+1];
int* mpot = new int[M.mu+1];
for(mrow=1;mrow<=M.mu;++mrow)
mnum[mrow] = 0;
for( t=1;t<=M.tu;++t)
++mnum[M.data[t].ii];
mpot[1] = 1;
//*(cpot+1)=1; //M中第一行第一個非零元在M中的序號為1
for(mrow=2;mrow<=M.mu;++mrow)
mpot[mrow] = mpot[mrow-1] + mnum[mrow-1];
//float si;
int po = 0;
float sr;
for (i=M.mu;i>=1;i--)
{
sr = r[i];
if(i<M.mu)
tp=mpot[i+1];
else
tp=M.tu+1;
for(int p=tp-1;p>mpot[i];p--)
{
t = M.data[p].jj;
sr = sr - M.data[p].e * solv[t];
}
po = mpot[i];
if (M.data[i].e==0)
{
solv[i] = 0;
}
else
solv[i] = sr/M.data[po].e;
}
}
void CMultiplyDlg::ICPCG()
{
float* ku = new float[K.nu+1];
float *z,*tz,*zp,*zpt;
float sum = 0;
int i,j;int cp = 0;
int count = 0;
z = new float[K.nu+1];
tz = new float[K.nu+1];
zp = new float[K.nu+1];
zpt= new float[K.nu+1];
float* kr = new float[K.mu+1];
float* krp = new float[K.mu+1];
float* krpt = new float[K.mu+1];
float rf,pt;
float *p = new float[K.nu+1];
K.vectMulSMat(U0,ku);
for (i=1;i<=K.mu;i++)
{
kr[i] = b[i] - ku[i];
sum+=kr[i]*kr[i];
}
TSMatrix ILt;
while (sum>0.0000001)
{
sum = 0;cp++;
offdiagsolve(ICL,kr,tz);
ICL.FastTrMatrix(ILt);
uppdiagsolve(ILt,tz,z);
count++;
if (count==1)
{
p = z;
zp = z;
krp = kr;
}
else
{
zpt = zp;
zp= z;
krpt = krp;
krp = kr;
float sup =0;
float supt = 0;
for (j=1;j<=K.mu;j++)
{
sup+= krp[j]*zp[j]; //[r(k-1)]T * z(k-1)
supt+=krpt[j]*zpt[j]; //[r(k-2)]T * z(k-2)
}
pt = sup/supt;
for (i=1;i<=K.mu;i++)
{
p[i] = zp[i] + pt*p[i];
}
}
float rsp = 0;
float pkp = 0;
int err = 0;
for (i=1;i<=K.mu;i++)
{
if (krp[i]<-1000 || zp[i]<-1000)
{
err = 1;
}
rsp+= krp[i]*zp[i]; //[r(k-1)]T * z(k-1)
}
float* kp = new float[K.mu+1];
K.vectMulSMat(p,kp); //kp = K*pk
for (i=1;i<=K.mu;i++)
{
pkp+=p[i]*kp[i]; //pkp = [p(k)]T *K * pk
}
rf = rsp/pkp;
for (i=1;i<=K.mu;i++)
{
kr[i] = krp[i] - rf*kp[i];
U0[i] = U0[i] + rf*p[i];
if (kr[i]!=0)
{
sum+=kr[i]*kr[i];
}
}
}
}
void CMultiplyDlg::iCholesky()
{
int *mnum = new int[K.mu+1];
int *mpot = new int[K.mu+1];
float* d = new float[K.mu+1];
int w = 3;
int h = 3;
indexSMat(K,mnum,mpot);
int i,j;
int tp;
int dk;
float sum = 0;
for(j=1;j<=K.mu;j++)
{
sum = 0;
tp = mpot[j];
//col = K.data[tp].jj;
for (i=tp;K.data[i].jj<j;i++)
{
dk = K.data[i].jj;
sum+= K.data[i].e*K.data[i].e/d[dk];
}
d[j] = K.data[i].e - sum;
}
TSMatrix dU,ofU,tsU;
dU.mu = K.mu;
dU.nu = K.nu;
dU.tu = K.mu;
dU.data = new Triple[dU.tu+1];
for (i=1;i<=K.mu;i++)
{
dU.data[i].e = d[i];
dU.data[i].ii = dU.data[i].jj = i;
}
ofU.data = new Triple[3+1];
ofU.mu = K.mu;
ofU.nu = K.nu;
ofU.tu = 3;
int gs = 0;
for (i=1;i<=K.mu;i++)
{
tp = mpot[i];
for (j=tp;K.data[j].jj<i;j++)
{
gs++;
ofU.data[gs] = K.data[j];
}
}
tsU.mu = dU.mu;
tsU.nu = dU.nu;
tsU.tu = dU.tu + ofU.tu;
tsU.data = new Triple[tsU.tu+1];
dU.addSMat(ofU,tsU);
int col;
for (i=1;i<=tsU.tu;i++)
{
col = tsU.data[i].jj;
tsU.data[i].e= 1/sqrtf(d[col])*tsU.data[i].e;
}
// ICL.mu = K.mu;
// ICL.nu = K.nu;
// ICL.tu = tsU.tu;
// ICL.data = new Triple[ICL.tu+1];
for (i=1;i<=ICL.tu;i++)
{
ICL.data[i] = tsU.data[i];
}
// delete tsU.data;
// delete dU.data;
// delete ofU.data;
// delete mnum;
// delete mpot;
// delete d;
// indexSMat(ICL,lmnum,lmpot);
}
void CMultiplyDlg::indexSMat(TSMatrix &M, int *mnum, int *mpot)
{
int mrow;
int t;
for(mrow=1;mrow<=M.mu;++mrow)
mnum[mrow] = 0;
for( t=1;t<=M.tu;++t)
++mnum[M.data[t].ii];
mpot[1] = 1;
//*(cpot+1)=1; //M中第一行第一個非零元在M中的序號為1
for(mrow=2;mrow<=M.mu;++mrow)
mpot[mrow] = mpot[mrow-1] + mnum[mrow-1];
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -