d9r5.cpp

《VC++ 常用数值算法》代码

#include "iostream.h"

void svdvar(double v[], int ma, int np, double w[], double cvm[], int ncvm)
{
	double wti[21];
	int i,j;
	for (i = 1; i<=ma; i++)
	{
		wti[i] = 0.0;
		if (w[i] != 0.0)
		{
			wti[i] = 1.0 / (w[i] * w[i]);
		}
	}
	for (i = 1; i<=ma; i++)
	{
		for (j = 1; j<=i; j++)
		{
			double sum1 = 0.0;
			for (int k = 1; k<=ma; k++)
			{
				sum1 = sum1 + v[(i-1)*np+k] * v[(j-1)*np+k] * wti[k];
			}
			cvm[(i-1)*ma+j] = sum1;
			cvm[(j-1)*ma+i] = sum1;
		}
	}
}

void main()
{
    //program d9r5
    //driver for routine svdvar
    int i,j,mp = 6;
    int ncvm,ma = 3;
    ncvm = ma;
    double v[37], w[7], cvm[10], tru[4][4];
    w[1] = 0.0;
    w[2] = 1.0;
    w[3] = 2.0;
    w[4] = 3.0;
    w[5] = 4.0;
    w[6] = 5.0;
    tru[1][1] = 1.25; tru[1][2] = 2.5; tru[1][3] = 3.75;
    tru[2][1] = 2.5;  tru[2][2] = 5.0; tru[2][3] = 7.5;
    tru[3][1] = 3.75; tru[3][2] = 7.5; tru[3][3] = 11.25;
    for (j = 1; j<=mp; j++)
	{
        for (i = 1; i<=mp; i++)
		{
            v[(i-1)*mp+j] = i;
		}
	}
    cout<<"Matrix v"<<endl;
    cout<<endl;
	cout.setf(ios::left);
    for (j = 1; j<=mp; j++)
	{
        for (i = 1; i<=mp; i++)
		{
			cout<<"   ";
			cout.width(6);
			cout<<v[(j-1)*mp+i];
        }
        cout<<endl;
    }
    cout<<endl;
    cout<<"Matrix w"<<endl;
    cout<<endl;
    for (i = 1; i<=mp; i++)
	{
		cout<<"   ";
		cout.width(6);
		cout<<w[i];
    }
    svdvar(v, ma, mp, w, cvm, ncvm);
    cout<<endl;cout<<endl;
    cout<<"Covariance matrix from svdvar"<<endl;
    cout<<endl;
    for (j = 1; j<=ma; j++)
	{
        for (i = 1; i<=ma; i++)
		{
			cout<<"   ";
			cout.width(12);
            cout<<cvm[(j-1)*ma+i];
        }
        cout<<endl;
    }
    cout<<endl;
    cout<<"Expected covariance matrix"<<endl;
    cout<<endl;
    for (j = 1; j<=ma; j++)
	{
        for (i = 1; i<=ma; i++)
		{
			cout<<"   ";
			cout.width(12);
            cout<<tru[j][i];
        }
        cout<<endl;
    }
}