multiplydlg.cpp

程序简单实现了对稀疏矩阵进行incomplete Cholesky分解的功能

// 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];
}