我真不认识小智
范文一：矩阵运算 数值转换
南京邮电大学通达学院程序设计(上机) 报告
题
目：
专业市场营销（物流管理）
学生姓名班级学号指导
教
师日
期
2010.6.14
评分项遵守机房规章制度
评
上机表现学习态度基础知识水平
分
程序设计能力课题功能实现情况算法设计合理性
细
用户界面友好性报告书写认真程度
则
报告内容完整性文字表达清楚程度问题回答准确性
评分等级备注
优秀良好中等及格不及格
教师签名：
年
月
日
注：评分等级分为五种：优秀、良好、中等、及格、不及
矩阵运算程序设计报告
一、问题描述
设计一个支持矩阵加减乘运算的程序
二课题分析
本题要求设计一个支持矩阵加减乘运算的程序，简单的功能框图如下：
矩阵运算
矩阵相加
矩阵相减矩阵相乘
输入矩阵输出矩阵
根据上图显示的功能要求，可以将问题转化为设计一个具有三个主要功能模块（矩阵相加、矩阵相减、矩阵相乘）的矩阵运算系统。同时，为了配合矩阵运算的需要，该矩阵运算系统必须包含输入和输出矩阵的功能。为简化运算量，这里，只考虑任意两个矩阵相加、相减和相乘的问题，三个或三个以上矩阵的运算并不执行。程序运行后，计划达到的目标如下：（1）提供可操作的主菜单：程序运行时，首先输出菜单，用于显示若干个可选
的功能选项。根据用户输入的选项来运行不同的功能，运行不同的函数。（2）输入信息：根据不同用户的需要，输入参与运算的矩阵，系统临时保存。（3）输出信息：将运算结果以矩阵的格式显现出来。
（4）矩阵相加：任意输入两个矩阵，如果满足相加条件，界面输出相加结果。
如果不满足相加条件，输出“矩阵不匹配”的提示，请用户重新操作。（5）矩阵相减：任意输入两个矩阵，如果满足相加条件，界面输出相减结果。
如果不满足相减条件，输出“矩阵不匹配”的提示，请用户重新操作。（6）矩阵相乘：任意输入两个矩阵，如果满足相乘条件，界面输出相乘结果。
如果不满足相乘条件，输出“矩阵不匹配”的提示，请用户重新操作。
三
{
数据结构的设计
class Str
private:
int m,n;//表示矩阵的行数和列数
double s[100][100];//双精度，定义二维数组存储矩阵public:int getm();
int getn();//调用私有数据void input();//输入矩阵函数void output();//输出矩阵函数
friend Str operator+(StrA,Str B); friend Str operator-(StrA,Str B);
friend Str operator*(StrA,Str B);//运算符的重载
};
矩阵类的各个主要部分具体设计方案如下：
（1）矩阵相加：首先判断要进行相加的两个矩阵是否匹配，即行数和列数是否相等，如果相等，则将两个矩阵类中对应位置的数组元素相加，结果保存到另一个矩阵类的数组中。
（2）矩阵相减：首先判断要进行相减的两个矩阵是否匹配，即行数和列数是否相等，如果相等，则将两个矩阵类中对应位置的数组元素相减，结果保存到另一个矩阵类的数组中。
（3）矩阵相乘：首先判断要进行相乘的两个矩阵是否匹配，即第一个矩阵的列数和第二个矩阵的行数是否相等，如果相等，根据矩阵相乘的定义，将两个矩阵类中对应位置的数组元素相乘后叠加，结果保存到另一个矩阵类的数组中。
为了使程序简洁明了，矩阵相加、相减和相乘函数均使用两个参数，用来表示参与运算的两个矩阵。同时，创建友元运算符函数operator()来实现“+”、“-”和“*”的重载，使其可以自由访问私有数据。
矩阵的匹配判断程序设计在主函数中，而不出现在矩阵相加、相减和相乘函数中。
主函数的设计：
（1）显示欢迎界面和选择菜单。
（2）使用switch()函数，使用户能够自由选择需要的运算。
（3）如果用户选择的选项不在菜单中，提示输入错误，并支持用户重新输入。
可以用default语句及循环实现。
（4）在执行某一选定的操作中，先调用类函数input（）来输入进行运算的矩
阵，然后进行矩阵是否匹配的判断。如果矩阵不匹配，马上终止运算并提示用户。如果矩阵匹配，则调用相应的重载运算符函数进行运算，并将运算结果用output（）函数输出。（5）利用while(i)不断循环，支持用户连续运算，仅当用户选择退出时，才
退出操作系统。
四处理结构的设计
（1）矩阵相加（减）函数
开始i=0
N
i
J=0
N
j
C.s[i][j]=A.s[i][j]+B.s[i][j]
j=j+1
C.m=A.m,C.n=A.n
对矩阵相减函数，此处换为相减后赋值。i=i+1
返回C
结束
矩阵相减函数的流程图与之类似。
(2)矩阵相乘函数
开始
i=0
N
i
j=0
N
j<>
C.s[i][j]=0,k=0N
j=j+1
k
Y
C.s[i][j]=C.s[i][j]+A.s[i][k]*B.s[k][j]
k=k+1
i=i+1
C.m=A.m,C.n=B.n
返回C
结束
五源程序
#include class Str { private: int m,n;//表示矩阵的行数和列数 double s[100][100];//定义二维数组存储矩阵public:int getm();int getn();void input();void output(); friend Str operator+(StrA,Str B); friend Str operator-(StrA,Str B); friend Str operator*(StrA,Str B);//运算符的重载}; Str operator+(StrA,Str B)//矩阵相加函数{ Str C; int i,j; for(i=0;i for(j=0;j C.s[i][j]=A.s[i][j]+B.s[i][j];C.m=A.m;C.n=A.n; return C; } Str operator-(StrA,Str B)//矩阵相减函数{ Str C; int i,j; for(i=0;i for(j=0;j C.s[i][j]=A.s[i][j]-B.s[i][j];C.m=A.m;C.n=A.n; return C; } Str operator*(StrA,Str B)//矩阵相乘函数 Str C; int i,j,k; for(i=0;i for(j=0;j<> C.s[i][j]=0; for(k=0;k C.s[i][j]=C.s[i][j]+A.s[i][k]*B.s[k][j];} C.m=A.m;C.n=B.n;return C; } void Str::input()//矩阵输入函数 { int i,j; do {cout<><> cin>>m>>n;}while(m<1||m>100||n<1||n>100);//保证输入矩阵的行数和 列数有效 cout<><><> for(j=0;j int i,j; cout<><><> for(j=0;j return m; int Str::getn(){ return n; } void menu(){ cout<><><><><><><><> int main(){ int choice,i=1;Str A,B,C; cout<><><><><><><><> while(i)//循环结构，使完成运算后自动进行下一次运算{ cout<> { case 1:cout<><> A.input(); cout<><> if(A.getm()==B.getm()&&A.getn()==B.getn())//判断是否符 合运算条件 { C=A+B;C.output();}else cout<> menu();//回到菜单，进行下一次运算 break;//使用break 语句，跳入下一个循环 case 2:cout<><> A.input(); cout<><> if(A.getm()==B.getm()&&A.getn()==B.getn()){ C=A-B;C.output();}else cout<<"矩阵不匹配！\n\n";menu();break; case=""><><> A.input(); cout<><> if(A.getn()==B.getm()){ C=A*B;C.output();}else cout<> case 4:cout<> cout<<"谢谢老师的辛勤教导我会继续努力！\n\n\n";i=0; 使其不符合循环条件，跳出循环break;=""><><> } }return 0; 六程序测试记录及调试记录 程序运行后，显示的主菜单界面： *******09005307宋丹丹设计******* 欢迎使用矩阵运算系统 1、矩阵相加运算 2、矩阵相减运算 3、矩阵相乘运算 4、退出运算系统 请选择： （1）矩阵相加测试： 输入：1 输出：矩阵A: 请输入矩阵的行数、列数： 输入：22 输出：请输入矩阵: 输入：1.12.2 3.34.4 输出：矩阵B: 请输入矩阵的行数、列数： 输入：22 输出：请输入矩阵: 输入：4.43.3 2.21.1 输出：得到的矩阵是： 5.55.5 5.55.5 （2）矩阵相减测试： 输入：2 输出：矩阵A: 请输入矩阵的行数、列数： 输入：22 输出：请输入矩阵: 输入：9.63.7 2.58.1 输出：矩阵B: 请输入矩阵的行数、列数： 输入：22 输出：请输入矩阵: 输入：3.31.8 0.73.5 输出：得到的矩阵是： 6.31.9 1.84.6 （3）矩阵相乘测试： 输入：3 输出：矩阵A: 请输入矩阵的行数、列数： 输入：22 输出：请输入矩阵: 输入：1.22 2.31 输出：矩阵B: 请输入矩阵的行数、列数： 输入：22 输出：请输入矩阵: 输入：3.23.1 2.34 输出：得到的矩阵是： 8.4411.72 9.6611.13 (4)退出矩阵运算系统测试 输入：4 输出：************谢谢使用！************ （5）在输入菜单选项时如果输入的内容不是1-4之间的数字，而是其他数字，系统将提示“Error Input! ”，若输入的为其他字符（字母或符号），也会出现这种结果。 （6）在输入的矩阵不匹配时，系统提示：矩阵不匹配！并重新显示菜单，请用户 重新选择。 七相关运行界面 （1）矩阵相加运算： （2）矩阵相减运算 （3）矩阵相乘运算 （4）矩阵不匹配测试 （5）退出矩阵运算系统 数制转换程序涉及报告 一问题描述 编制不同数制间的转换程序。要求提供输入输出界面，当输完一个任意十进制数字，程序能自动将其转化为另外的数制表示的值，包括二进制、八进制和十六进制，其中转化用算法实现，而不是用printf 函数显示。 二课题分析 将十进制数转换为其他进制的数，其实质是以将被换为的进制数作为除数，从给出的数起依次作除求余，直至整除后结果为零，然后将运算中所得的余数按倒序排列即得到数制转换后所得 的数。 在定义数组时需要注意到无符号长整型数的取值范围，提供足够但有效地数组长度。 转换为二进制的数与转换为八进制的数的算法基本一致，即依次记下求余所得，并判断整除后是否为零决定是否继续运算，而转换为十六进制的数时需要注意前十个为数字表示，而后面六个数为用大写的’A’~’F’表示。 运算结束后需要按照从后向前的顺序依次输出，此时需要注意因转换运算中的’i++’，此时数组最后一个位置没有数据存储，需要回到前一个位置。 考虑到充分满足使用者的需求，完成一次运算后自动询问是否继续，继续则继续运算，否则即退出系统。 三数据结构的设计 unsigned long a,temp;//存储所要转换的数和整除后的数char b[64];//存储二进制数的数组 char o[21];//存储十进制数的数组 char h[16];//存储十六进制数的数组 在此，所定义的三个数组的长度即分别有效地考虑到无符号长整型数的取值范围 b[i]=temp%2+'0';//求余；强制类型转换为实数 temp/=2;//整除运算 在此需根据temp 的值判断是否运算完毕 while(i>=0) {printf("%c",b[i]); i--;//从后向前依次输出 } while(temp&&i<> { int x=temp%16; if(x<> else h[i]=x-10+'A';//前十个为数字，后六个为字母temp/=16; i++; } int m=1; while(m) {…… scanf("%d",&m); switch(m) {case1:continue; case 0:printf("************谢谢使用本系统 ************"); }} While 循环的设计可有效提高系统的方便性和实用性。四处理结构的设计 主函数 开始 输入待转换数进行数制转换运算输出转换后的数 结束 二进制转换部分函数 开始 temp=a;i=0b[i]=temp%2+’0’temp/=2;i++temp!=0 &&i<> ？ 否是 结束 八进制与之类似 十六进制转换部分函数 开始 temp=a;i=0 i<> 否 是 b[i]=temp%2+’0’b[i]=temp%2+’0’ temp/=2;i++ temp!=0 &&i<> ? 是 否 结束 输出部分函数 开始 输出b[i] i-- i>=0 否是 结束 五源程序 #include int main() { unsigned long a,temp;//定义无符号长整型 char b[64]; char o[21]; char h[16];//无符号长整型数的取值范围 int i,m=1; printf("*********09005307宋丹丹设计*********\n\n\n"); while(m) { printf("Inputnumber:"); scanf("%ld",&a); temp=a; i=0; while(temp&&i<> { b[i]=temp%2+'0';//求余；强制类型转换为实数 temp/=2;//整除运算 i++; } printf("\n二进制:"); i--;//返回最后一个数的位置 while(i>=0) { printf("%c",b[i]); i--;//从后向前依次输出 } printf("\n"); temp=a; i=0; while(temp&&i<> { o[i]=temp%8+'0'; temp/=8; i++; } printf("八进制: i--; while(i>=0) { printf("%c",o[i]); i--; } printf("\n");"); temp=a; i=0; while(temp&&i<> { int x=temp%16; if(x<> else h[i]=x-10+'A';//前十个为数字，后六个为字母temp/=16; i++; } printf("十六进制:"); i--; while(i>=0&&i<> { printf("%c",h[i]); i--; } printf("\n\n\n"); do{ printf("1、继续\n0、退出\n请选择：\n"); scanf("%d",&m); }while(m<0||m>1); switch(m) { case 1:continue; case 0:printf("************谢谢使用 ************"); }本系统 } printf("\n\n\n"); printf("**谢谢老师的辛勤教导我会继续努力！**\n\n\n");return 0; } 六调试记录 程序运行后，输出界面： *********09005307宋丹丹设计********* Input number: 第一次输入：65 输出： 二进制: 八进制:1000001101 十六进制:41 1、继续 0、退出 输入：1 输出：Inputnumber: 第二次输入：165 输出： 二进制: 八进制:10100101245 十六进制:A5 1、继续 0、退出 输入：0 输出：************谢谢使用本系统************七 相关界面 卷积的数值运算 1、主函数（计算任意两个函数的卷积） %%%用数值计算方法计算连续信号卷积积分f(t)=f1(t)*f2(t)%%% %k是指f(t)对应的时间向量；p 是指取样时间间隔； %f1是指f1(t)的非零样值向量；f2是指f2(t)的非零样值向量； %k1是指f1(t)的对应时间向量；k2是指f2(t)的对应时间向量； function[f,k]=sconv(f1,f2,k1,k2,p) f=conv(f1,f2); f=f*p; k0=k1(1)+k2(1); %k0为f(t)时间向量的起点位置 k3=length(f1)+length(f2)-2; %k3为f(t)的样值总宽度 k=k0:p:k3*p; subplot(2,2,1);plot(k1,f1); title('f1(t)');xlabel('t');ylabel('f(1)'); %绘制f1(t)的时域图 subplot(2,2,2);plot(k2,f2); title('f2(t)');xlabel('t');ylabel('f(2)'); %绘制f2(t)的时域图 subplot(2,2,3);plot(k,f); %绘制f(t)的时域波形图 h=get(gca,'position');h(3)=2*h(3); set(gca,'position',h); %将f(t)函数的横坐标范围扩为原图的2倍 title('f(t)=f1(t)*f2(t)'); xlabel('t'); ylabel('f(t)'); 2、子函数（计算f1=cos(k1)和f2=sin(k1)两个函数的卷积）： p=0.1; k1=0:p:20; f1=cos(k1); k2=k1; f2=sin(k1); [f,k]=sconv(f1,f2,k1,k2,p);%调用主函数 3、MATLAB 实现效果： 第 1 页 数值计算程序大放送-矩阵运算 ////////////////////////////////////////////////////////////// //实矩阵相乘 //计算矩阵A(m*n)和B(n*k)的乘积，结果保存在C(m*k)中 //a-长度为m*n的数组 //b-长度为n*k的数组 //c-长度为m*k的数组，存放结果 void damul(double a[],double b[],int m,int n,int k,double c[]); ////////////////////////////////////////////////////////////// //计算矩阵A(m*n)的转置矩阵AT(n*m)和B(m*k)的乘积，结果保存在C(n*k)中 //添加的函数，非原书程序 //a-长度为m*n的数组 //b-长度为m*k的数组 //c-长度为n*k的数组，存放结果 void ATdotB(double a[],double b[],int m,int n,int k,double c[]); ////////////////////////////////////////////////////////////// //计算矩阵A(m*n)和B(k*n)的转置矩阵BT(n*k)的乘积，结果保存在C(m*k)中 //添加的函数，非原书程序 //a-长度为m*n的数组 //b-长度为k*n的数组 //c-长度为m*k的数组，存放结果 void AdotBT(double a[],double b[],int m,int n,int k,double c[]); ////////////////////////////////////////////////////////////// //实矩阵求逆 //全选主元高斯-约当法 //a-长度为n*n的数组, n*n矩阵 //n 矩阵的维数 int dcinv(double a[],int n); ////////////////////////////////////////////////////////////// //对称正定矩阵求逆 //a-长度为n*n的数组, n*n矩阵 //n 矩阵的维数 int desgj(double a[],int n); ////////////////////////////////////////////////////////////// //托伯利兹(Toeplitz)矩阵求逆的特兰持(Trench)方法 //t-长度为n的数组,存放n阶T型矩阵中的上三角元素t0,t1,t2...tn-1 //tt-长度为n的数组,从tt[1]开始依次存放tt[1]...tt[n-1] //n-矩阵的阶数 //b-长度为n*n的数组，返回时存放逆矩阵 int dftrn(double t[],double tt[],int n,double b[]); ////////////////////////////////////////////////////////////// //求矩阵的行列式值 //全选主元高斯消去法 //a-长度为n*n的数组 //n-矩阵的阶数 double dhdet(double a[],int n); ////////////////////////////////////////////////////////////// //对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值 //返回值小于0表示程序工作失败，还输出"fail"; //返回值大于0表示正常返回 //a-长度为n*n的数组，存放正定矩阵， // 返回时下三角部分存放分解后的下三角矩阵L,其余元素为0 //n-正定矩阵的阶数 //det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值 int dicll(double a[],int n,double *det); ////////////////////////////////////////////////////////////// //矩阵的三角分解(LU) //其中下三角阵L的主对角元素为1。 //a-长度为N*N的矩阵，返回时为L+U-I //n-矩阵的阶数 //l-返回下三角矩阵 //u-返回上三角矩阵 int djlu(double a[],int n,double l[],double u[]); ////////////////////////////////////////////////////////////// //实数矩阵的QR分解法 //用Householder变换对一般m*n阶的实数矩阵进行QR分解 //a-长度为m*n的一维数组,返回时其左上三角部分存放上三角矩阵R //m-矩阵的行数 //n-矩阵的列数 //q-长度为m*m的矩阵，返回时存放正交矩阵Q int dkqr(double a[],int m,int n,double q[]); ////////////////////////////////////////////////////////////// //奇异值分解法求广义逆 //本函数返回值小于0表示在奇异值分解过程, //中迭代值超过了60次还未满足精度要求. //返回值大于0表示正常返回。 //a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0 //m-矩阵的行数 //n-矩阵的列数 //aa-长度为n*m的数组，返回式存放A的广义逆 //eps-精度要求 //u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U //v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V //ka-整型变量，其值为max(n,m)+1 //调用函数:dluav() int dginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka); ////////////////////////////////////////////////////////////// //实数矩阵的奇异值分解 //利用Householder变换及变形QR算法 //a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0 //m-矩阵的行数 //n-矩阵的列数 //u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U //v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V //eps-精度要求 //ka-整型变量，其值为max(n,m)+1 //调用函数:dluav(),ppp(),sss() int dluav(double a[],int m,int n,double u[],double v[],double eps,int ka); ////////////////////////////////////////////////////////////// #include "math.h" #include "stdio.h" //实矩阵相乘 //计算矩阵A(m*n)和B(n*k)的乘积，结果保存在C(m*k)中 //a-长度为m*n的数组 //b-长度为n*k的数组 //c-长度为m*k的数组，存放结果 void damul(double a[],double b[],int m,int n,int k,double c[]) { int i,j,l,u; for (i=0; i for (j=0; j { u=i*k+j; c[u]=0.0; for (l=0; l { c[u]+=a[i*n+l]*b[l*k+j]; } } } return; } //计算矩阵A(m*n)的转置矩阵AT(n*m)和B(m*k)的乘积，结果保存在C(n*k)中 //添加的函数，非原书程序 //a-长度为m*n的数组 //b-长度为m*k的数组 //c-长度为n*k的数组，存放结果 void ATdotB(double a[],double b[],int m,int n,int k,double c[]) { int i,j,l,u; for (i=0; i { for (j=0; j { u=i*k+j; c[u]=0.0; for (l=0; l { c[u]+=a[l*n+i]*b[l*k+j]; } } } return; } //计算矩阵A(m*n)和B(k*n)的转置矩阵BT(n*k)的乘积，结果保存在C(m*k)中 //添加的函数，非原书程序 //a-长度为m*n的数组 //b-长度为k*n的数组 //c-长度为m*k的数组，存放结果 void AdotBT(double a[],double b[],int m,int n,int k,double c[]) { int i,j,l,u; for (i=0; i for (j=0; j { u=i*k+j; c[u]=0.0; for (l=0; l { c[u]+=a[i*n+l]*b[j*n+l]; } } } return; } //实矩阵求逆 //a-长度为n*n的数组, n*n矩阵 //n 矩阵的维数 int dcinv(double a[],int n) { int *is,*js,i,j,k,l,u,v; double d,p; is=new int[n]; is=malloc(n*sizeof(int)); js=malloc(n*sizeof(int)); for (k=0; k<=n-1; k++)="" {=""> d=0.0; for (i=k; i<=n-1; i++)=""> { for (j=k; j<=n-1; j++)=""> { l=i*n+j; p=fabs(a[l]); if (p>d) { d=p; is[k]=i; js[k]=j; } } } if (d+1.0==1.0) { free(is); free(js); printf("err**not inv\n"); return(0); } if (is[k]!=k) { for (j=0; j<=n-1; j++)=""> { u=k*n+j; v=is[k]*n+j; p=a[u]; a[u]=a[v]; a[v]=p; } } if (js[k]!=k) { for (i=0; i<=n-1; i++)=""> { u=i*n+k; v=i*n+js[k]; p=a[u]; a[u]=a[v]; a[v]=p; } } l=k*n+k; a[l]=1.0/a[l]; for (j=0; j<=n-1; j++)=""> { if (j!=k) { u=k*n+j; a[u]=a[u]*a[l]; } } for (i=0; i<=n-1; i++)=""> { if (i!=k) { for (j=0; j<=n-1; j++)=""> { if (j!=k) { u=i*n+j; a[u]=a[u]-a[i*n+k]*a[k*n+j]; } } } } for (i=0; i<=n-1; i++)=""> { if (i!=k) { u=i*n+k; a[u]=-a[u]*a[l]; } } } for (k=n-1; k>=0; k--) { if (js[k]!=k) { for (j=0; j<=n-1; j++)=""> { u=k*n+j; v=js[k]*n+j; p=a[u]; a[u]=a[v]; a[v]=p; } } if (is[k]!=k) { for (i=0; i<=n-1; i++)=""> { u=i*n+k; v=i*n+is[k]; p=a[u]; a[u]=a[v]; a[v]=p; } } } free(is); free(js); return(1); //对称正定矩阵求逆 //a-长度为n*n的数组, n*n矩阵 //n 矩阵的维数 int desgj(double a[],int n) { int i,j,k,m; double w,g,*b; b=malloc(n*sizeof(double)); for (k=0; k<=n-1; k++)="" {=""> w=a[0]; if (fabs(w)+1.0==1.0) { free(b); printf("fail\n"); return(-2); } m=n-k-1; for (i=1; i<=n-1; i++)=""> { g=a[i*n]; b[i]=g/w; if (i<=m)> { b[i]=-b[i]; } for (j=1; j<=i; j++)=""> { a[(i-1)*n+j-1]=a[i*n+j]+g*b[j]; } } a[n*n-1]=1.0/w; for (i=1; i<=n-1; i++)=""> { a[(n-1)*n+i-1]=b[i]; } } for (i=0; i<=n-2; i++)="" {=""> for (j=i+1; j<=n-1; j++)=""> { a[i*n+j]=a[j*n+i]; } } free(b); return(2); } //托伯利兹(Toeplitz)矩阵求逆的特兰持(Trench)方法 //t-长度为n的数组,存放n阶T型矩阵中的上三角元素t0,t1,t2...tn-1 //tt-长度为n的数组,从tt[1]开始依次存放tt[1]...tt[n-1] //n-矩阵的阶数 //b-长度为n*n的数组，返回时存放逆矩阵 int dftrn(double t[],double tt[],int n,double b[]) { int i,j,k; double a,s,*c,*r,*p; c=malloc(n*sizeof(double)); r=malloc(n*sizeof(double)); p=malloc(n*sizeof(double)); if (fabs(t[0])+1.0==1.0) { free(c); free(r); free(p); printf("fail\n"); return(-1); } a=t[0]; c[0]=tt[1]/t[0]; r[0]=t[1]/t[0]; for (k=0; k<=n-3; k++)="" {=""> s=0.0; for (j=1; j<=k+1; j++)=""> { s=s+c[k+1-j]*tt[j]; } s=(s-tt[k+2])/a; for (i=0; i<=k;> { p[i]=c[i]+s*r[k-i]; } c[k+1]=-s; s=0.0; for (j=1; j<=k+1; j++)=""> { s=s+r[k+1-j]*t[j]; } s=(s-t[k+2])/a; for (i=0; i<=k; i++)=""> { r[i]=r[i]+s*c[k-i]; c[k-i]=p[k-i]; } r[k+1]=-s; a=0.0; for (j=1; j<=k+2; j++)=""> { a=a+t[j]*c[j-1]; } a=t[0]-a; if (fabs(a)+1.0==1.0) { free(c); free(r); free(p); printf("fail\n"); return(-1); } } b[0]=1.0/a; for (i=0; i<=n-2; i++)="" {=""> k=i+1; j=(i+1)*n; b[k]=-r[i]/a; b[j]=-c[i]/a; } for (i=0; i<=n-1; i++)="" {=""> for (j=0; j<=n-2; j++)=""> { k=(i+1)*n+j+1; b[k]=b[i*n+j]-c[i]*b[j+1]; b[k]=b[k]+c[n-j-2]*b[n-i-1]; } } free(c); free(r); free(p); return(1); } //求矩阵的行列式值 //全选主元高斯消去法 //a-长度为n*n的数组 //n-矩阵的阶数 double dhdet(double a[],int n) { int i,j,k,is,js,l,u,v; double f,det,q,d; f=1.0; det=1.0; for (k=0; k<=n-2; k++)="" {=""> q=0.0; for (i=k; i<=n-1; i++)=""> for (j=k; j<=n-1; j++)=""> { l=i*n+j; d=fabs(a[l]); if (d>q) { q=d; is=i; js=j; } } if (q+1.0==1.0) { det=0.0; return(det); } if (is!=k) { f=-f; for (j=k; j<=n-1; j++)=""> { u=k*n+j; v=is*n+j; d=a[u]; a[u]=a[v]; a[v]=d; } } if (js!=k) { f=-f; for (i=k; i<=n-1; i++)=""> { u=i*n+js; v=i*n+k; d=a[u]; a[u]=a[v]; a[v]=d; } } l=k*n+k; det=det*a[l]; for (i=k+1; i<=n-1; i++)=""> { d=a[i*n+k]/a[l]; for (j=k+1; j<=n-1; j++)=""> { u=i*n+j; a[u]=a[u]-d*a[k*n+j]; } } } det=f*det*a[n*n-1]; return(det); } //对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值 //返回值小于0表示程序工作失败，还输出"fail"; //返回值大于0表示正常返回 //a-长度为n*n的数组，存放正定矩阵， // 返回时下三角部分存放分解后的下三角矩阵L,其余元素为0 //n-正定矩阵的阶数 //det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值 int dicll(double a[],int n,double *det) { int i,j,k,u,v,l; double d; if ((a[0]+1.0==1.0)(a[0]<0.0))> { printf("fail\n"); return(-2); } a[0]=sqrt(a[0]); d=a[0]; for (i=1; i<=n-1; i++)="" {=""> u=i*n; a[u]=a[u]/a[0]; } for (j=1; j<=n-1; j++)="" {=""> l=j*n+j; for (k=0; k<=j-1; k++)=""> { u=j*n+k; a[l]=a[l]-a[u]*a[u]; } if ((a[l]+1.0==1.0)(a[l]<0.0))> { printf("fail\n"); return(-2); } a[l]=sqrt(a[l]); d=d*a[l]; for (i=j+1; i<=n-1; i++)=""> { u=i*n+j; for (k=0; k<=j-1; k++)=""> { a[u]=a[u]-a[i*n+k]*a[j*n+k]; } a[u]=a[u]/a[l]; } } *det=d*d; for (i=0; i<=n-2; i++)="" {=""> for (j=i+1; j<=n-1; j++)=""> { a[i*n+j]=0.0; } } return(2); } //求矩阵的行列式值 //全选主元高斯消去法 //a-长度为n*n的数组 //n-矩阵的阶数 double dhdet(double a[],int n) { int i,j,k,is,js,l,u,v; double f,det,q,d; f=1.0; det=1.0; for (k=0; k<=n-2; k++)="" {=""> q=0.0; for (i=k; i<=n-1; i++)=""> for (j=k; j<=n-1; j++)=""> { l=i*n+j; d=fabs(a[l]); if (d>q) { q=d; is=i; js=j; } } if (q+1.0==1.0) { det=0.0; return(det); } if (is!=k) { f=-f; for (j=k; j<=n-1; j++)=""> { u=k*n+j; v=is*n+j; d=a[u]; a[u]=a[v]; a[v]=d; } } if (js!=k) { f=-f; for (i=k; i<=n-1; i++)=""> { u=i*n+js; v=i*n+k; d=a[u]; a[u]=a[v]; a[v]=d; } } l=k*n+k; det=det*a[l]; for (i=k+1; i<=n-1; i++)=""> { d=a[i*n+k]/a[l]; for (j=k+1; j<=n-1; j++)=""> { u=i*n+j; a[u]=a[u]-d*a[k*n+j]; } } } det=f*det*a[n*n-1]; return(det); } //对称正定矩阵的乔里斯基(Cholesky)分解与行列式求值 //返回值小于0表示程序工作失败，还输出"fail"; //返回值大于0表示正常返回 //a-长度为n*n的数组，存放正定矩阵， // 返回时下三角部分存放分解后的下三角矩阵L,其余元素为0 //n-正定矩阵的阶数 //det-指向双精度实型变量的指针，返回时该指针指向的变量存放行列式的值 int dicll(double a[],int n,double *det) { int i,j,k,u,v,l; double d; if ((a[0]+1.0==1.0)(a[0]<0.0))> { printf("fail\n"); return(-2); } a[0]=sqrt(a[0]); d=a[0]; for (i=1; i<=n-1; i++)="" {=""> u=i*n; a[u]=a[u]/a[0]; } for (j=1; j<=n-1; j++)="" {=""> l=j*n+j; for (k=0; k<=j-1;> { u=j*n+k; a[l]=a[l]-a[u]*a[u]; } if ((a[l]+1.0==1.0)(a[l]<0.0))> { printf("fail\n"); return(-2); } a[l]=sqrt(a[l]); d=d*a[l]; for (i=j+1; i<=n-1; i++)=""> { u=i*n+j; for (k=0; k<=j-1; k++)=""> { a[u]=a[u]-a[i*n+k]*a[j*n+k]; } a[u]=a[u]/a[l]; } } *det=d*d; for (i=0; i<=n-2; i++)="" {=""> for (j=i+1; j<=n-1; j++)=""> { a[i*n+j]=0.0; } } return(2); } //矩阵的三角分解(LU) //其中下三角阵L的主对角元素为1。 //a-长度为N*N的矩阵，返回时为L+U-I //n-矩阵的阶数 //l-返回下三角矩阵 //u-返回上三角矩阵 int djlu(double a[],int n,double l[],double u[]) { int i,j,k,w,v,ll; for (k=0; k<=n-2;> { ll=k*n+k; if (fabs(a[ll])+1.0==1.0) { printf("fail\n"); return(0); } for (i=k+1; i<=n-1; i++)=""> { w=i*n+k; a[w]=a[w]/a[ll]; } for (i=k+1; i<=n-1; i++)=""> { w=i*n+k; for (j=k+1; j<=n-1; j++)=""> { v=i*n+j; a[v]=a[v]-a[w]*a[k*n+j]; } } } for (i=0; i<=n-1; i++)="" {=""> for (j=0; j { w=i*n+j; l[w]=a[w]; u[w]=0.0; } w=i*n+i; l[w]=1.0; u[w]=a[w]; for (j=i+1; j<=n-1; j++)=""> { w=i*n+j; l[w]=0.0; u[w]=a[w]; } } return(1); } //实数矩阵的QR分解法 //用Householder变换对一般m*n阶的实数矩阵进行QR分解 //a-长度为m*n的一维数组,返回时其左上三角部分存放上三角矩阵R //m-矩阵的行数 //n-矩阵的列数 //q-长度为m*m的矩阵，返回时存放正交矩阵Q int dkqr(double a[],int m,int n,double q[]) { int i,j,k,l,nn,p,jj; double u,alpha,w,t; if (m { printf("fail\n"); return(0); } for (i=0; i<=m-1; i++)=""> { for (j=0; j<=m-1; j++)=""> { l=i*m+j; q[l]=0.0; if (i==j) { q[l]=1.0; } } } nn=n; if (m==n) { nn=m-1; } for (k=0; k<=nn-1; k++)="" {=""> u=0.0; l=k*n+k; for (i=k; i<=m-1; i++)=""> { w=fabs(a[i*n+k]); if (w>u) { u=w; } } alpha=0.0; for (i=k; i<=m-1; i++)=""> { t=a[i*n+k]/u; alpha=alpha+t*t; } if (a[l]>0.0) { u=-u; } alpha=u*sqrt(alpha); if (fabs(alpha)+1.0==1.0) { printf("fail\n"); return(0); } u=sqrt(2.0*alpha*(alpha-a[l])); if ((u+1.0)!=1.0) { a[l]=(a[l]-alpha)/u; for (i=k+1; i<=m-1; i++)=""> { p=i*n+k; a[p]=a[p]/u; } for (j=0; j<=m-1; j++)=""> { t=0.0; for (jj=k; jj<=m-1; jj++)=""> { t=t+a[jj*n+k]*q[jj*m+j]; } for (i=k; i<=m-1; i++)=""> { p=i*m+j; q[p]=q[p]-2.0*t*a[i*n+k]; } } for (j=k+1; j<=n-1; j++)=""> { t=0.0; for (jj=k; jj<=m-1; jj++)=""> { t=t+a[jj*n+k]*a[jj*n+j]; } for (i=k; i<=m-1; i++)=""> { p=i*n+j; a[p]=a[p]-2.0*t*a[i*n+k]; } } a[l]=alpha; for (i=k+1; i<=m-1; i++)=""> { a[i*n+k]=0.0; } } } for (i=0; i<=m-2; i++)="" {=""> for (j=i+1; j<=m-1;j++)> { p=i*m+j; l=j*m+i; t=q[p]; q[p]=q[l]; q[l]=t; } } return(1); } //奇异值分解法求广义逆 //本函数返回值小于0表示在奇异值分解过程, //中迭代值超过了60次还未满足精度要求. //返回值大于0表示正常返回。 //a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0 //m-矩阵的行数 //n-矩阵的列数 //aa-长度为n*m的数组，返回式存放A的广义逆 //eps-精度要求 //u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U //v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V //ka-整型变量，其值为max(n,m)+1 //调用函数:dluav() int dginv(double a[],int m,int n,double aa[],double eps,double u[],double v[],int ka) { int i,j,k,l,t,p,q,f; i=dluav(a,m,n,u,v,eps,ka); if (i<0)> { return(-1); } j=n; if (m { j=m; } j=j-1; k=0; while ((k<=j)&&(a[k*n+k]!=0.0))> { k=k+1; } k=k-1; for (i=0; i<=n-1; i++)="" {=""> for (j=0; j<=m-1; j++)=""> { t=i*m+j; aa[t]=0.0; for (l=0; l<=k; l++)=""> { f=l*n+i; p=j*m+l; q=l*n+l; aa[t]=aa[t]+v[f]*u[p]/a[q]; } } } return(1); } //实数矩阵的奇异值分解 //利用Householder变换及变形QR算法 //a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0 //m-矩阵的行数 //n-矩阵的列数 //u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U //v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V //eps-精度要求 //ka-整型变量，其值为max(n,m)+1 //调用函数:dluav(),ppp(),sss() static void ppp(double a[],double e[],double s[],double v[],int m,int n); static void sss(double fg[2],double cs[2]); int dluav(double a[],int m,int n,double u[],double v[],double eps,int ka) { int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks; double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2]; double *s,*e,*w; s=malloc(ka*sizeof(double)); e=malloc(ka*sizeof(double)); w=malloc(ka*sizeof(double)); it=60; k=n; if (m-1 { k=m-1; } l=m; if (n-2 { l=n-2; } if (l<0)> { l=0; } ll=k; if (l>k) { ll=l; } if (ll>=1) { for (kk=1; kk<=ll; kk++)=""> { if (kk<=k)> { d=0.0; for (i=kk; i<=m; i++)=""> { ix=(i-1)*n+kk-1; d=d+a[ix]*a[ix]; } s[kk-1]=sqrt(d); if (s[kk-1]!=0.0) { ix=(kk-1)*n+kk-1; if (a[ix]!=0.0) { s[kk-1]=fabs(s[kk-1]); if (a[ix]<0.0)> { s[kk-1]=-s[kk-1]; } } for (i=kk; i<=m; i++)=""> { iy=(i-1)*n+kk-1; a[iy]=a[iy]/s[kk-1]; } a[ix]=1.0+a[ix]; } s[kk-1]=-s[kk-1]; } if (n>=kk+1) { for (j=kk+1; j<=n; j++)=""> { if ((kk<=k)&&(s[kk-1]!=0.0))> { d=0.0; for (i=kk; i<=m; i++)=""> { ix=(i-1)*n+kk-1; iy=(i-1)*n+j-1; d=d+a[ix]*a[iy]; } d=-d/a[(kk-1)*n+kk-1]; for (i=kk; i<=m; i++)=""> { ix=(i-1)*n+j-1; iy=(i-1)*n+kk-1; a[ix]=a[ix]+d*a[iy]; } } e[j-1]=a[(kk-1)*n+j-1]; } } if (kk<> { for (i=kk; i<=m; i++)=""> { ix=(i-1)*m+kk-1; iy=(i-1)*n+kk-1; u[ix]=a[iy]; } } if (kk<=l)> { d=0.0; for (i=kk+1; i<=n; i++)=""> { d=d+e[i-1]*e[i-1]; } e[kk-1]=sqrt(d); if (e[kk-1]!=0.0) { if (e[kk]!=0.0) { e[kk-1]=fabs(e[kk-1]); if (e[kk]<0.0)> { e[kk-1]=-e[kk-1]; } } for (i=kk+1; i<=n; i++)=""> { e[i-1]=e[i-1]/e[kk-1]; } e[kk]=1.0+e[kk]; } e[kk-1]=-e[kk-1]; if ((kk+1<=m)&&(e[kk-1]!=0.0))> { for (i=kk+1; i<=m; i++)=""> { w[i-1]=0.0; } for (j=kk+1; j<=n; j++)=""> { for (i=kk+1; i<=m; i++)=""> { w[i-1]=w[i-1]+e[j-1]*a[(i-1)*n+j-1]; } } for (j=kk+1; j<=n; j++)=""> { for (i=kk+1; i<=m; i++)=""> { ix=(i-1)*n+j-1; a[ix]=a[ix]-w[i-1]*e[j-1]/e[kk]; } } } for (i=kk+1; i<=n; i++)=""> { v[(i-1)*n+kk-1]=e[i-1]; } } } } mm=n; if (m+1 { mm=m+1; } if (k { s[k]=a[k*n+k]; } if (m { s[mm-1]=0.0; } if (l+1 { e[l]=a[l*n+mm-1]; } e[mm-1]=0.0; nn=m; if (m>n) { nn=n; } if (nn>=k+1) { for (j=k+1; j<=nn; j++)=""> { for (i=1; i<=m; i++)=""> { u[(i-1)*m+j-1]=0.0; } u[(j-1)*m+j-1]=1.0; } } if (k>=1) { for (ll=1; ll<=k; ll++)=""> { kk=k-ll+1; iz=(kk-1)*m+kk-1; if (s[kk-1]!=0.0) { if (nn>=kk+1) { for (j=kk+1; j<=nn; j++)=""> { d=0.0; for (i=kk; i<=m; i++)="" {=""> ix=(i-1)*m+kk-1; iy=(i-1)*m+j-1; d=d+u[ix]*u[iy]/u[iz]; } d=-d; for (i=kk; i<=m; i++)=""> { ix=(i-1)*m+j-1; iy=(i-1)*m+kk-1; u[ix]=u[ix]+d*u[iy]; } } } for (i=kk; i<=m; i++)=""> { ix=(i-1)*m+kk-1; u[ix]=-u[ix]; } u[iz]=1.0+u[iz]; if (kk-1>=1) { for (i=1; i<=kk-1; i++)=""> { u[(i-1)*m+kk-1]=0.0; } } } else { for (i=1; i<=m; i++)=""> { u[(i-1)*m+kk-1]=0.0; } u[(kk-1)*m+kk-1]=1.0; } } } for (ll=1; ll<=n; ll++)="" {=""> kk=n-ll+1; iz=kk*n+kk-1; if ((kk<=l)&&(e[kk-1]!=0.0))> { for (j=kk+1; j<=n; j++)=""> { d=0.0; for (i=kk+1; i<=n; i++)="" {=""> ix=(i-1)*n+kk-1; iy=(i-1)*n+j-1; d=d+v[ix]*v[iy]/v[iz]; } d=-d; for (i=kk+1; i<=n; i++)=""> { ix=(i-1)*n+j-1; iy=(i-1)*n+kk-1; v[ix]=v[ix]+d*v[iy]; } } } for (i=1; i<=n; i++)=""> { v[(i-1)*n+kk-1]=0.0; } v[iz-n]=1.0; } for (i=1; i<=m; i++)="" {=""> for (j=1; j<=n; j++)=""> { a[(i-1)*n+j-1]=0.0; } } m1=mm; it=60; while (1==1) { if (mm==0) { ppp(a,e,s,v,m,n); free(s); free(e); free(w); return(1); } if (it==0) { ppp(a,e,s,v,m,n); free(s); free(e); free(w); return(-1); } kk=mm-1; while ((kk!=0)&&(fabs(e[kk-1])!=0.0)) { d=fabs(s[kk-1])+fabs(s[kk]); dd=fabs(e[kk-1]); if (dd>eps*d) { kk=kk-1; } else { e[kk-1]=0.0; } } if (kk==mm-1) { kk=kk+1; if (s[kk-1]<0.0)> { s[kk-1]=-s[kk-1]; for (i=1; i<=n; i++)=""> { ix=(i-1)*n+kk-1; v[ix]=-v[ix]; } } while ((kk!=m1)&&(s[kk-1] { d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d; if (kk { for (i=1; i<=n; i++)=""> { ix=(i-1)*n+kk-1; iy=(i-1)*n+kk; d=v[ix]; v[ix]=v[iy]; v[iy]=d; } } if (kk { for (i=1; i<=m; i++)=""> { ix=(i-1)*m+kk-1; iy=(i-1)*m+kk; d=u[ix]; u[ix]=u[iy]; u[iy]=d; } } kk=kk+1; } it=60; mm=mm-1; } else { ks=mm; while ((ks>kk)&&(fabs(s[ks-1])!=0.0)) { d=0.0; if (ks!=mm) { d=d+fabs(e[ks-1]); } if (ks!=kk+1) { d=d+fabs(e[ks-2]); } dd=fabs(s[ks-1]); if (dd>eps*d) { ks=ks-1; } else { s[ks-1]=0.0; } } if (ks==kk) { kk=kk+1; d=fabs(s[mm-1]); t=fabs(s[mm-2]); if (t>d) { d=t; } t=fabs(e[mm-2]); if (t>d) { d=t; } t=fabs(s[kk-1]); if (t>d) { d=t; } t=fabs(e[kk-1]); if (t>d) { d=t; } sm=s[mm-1]/d; sm1=s[mm-2]/d; em1=e[mm-2]/d; sk=s[kk-1]/d; ek=e[kk-1]/d; b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0; c=sm*em1; c=c*c; shh=0.0; if ((b!=0.0)(c!=0.0)) { shh=sqrt(b*b+c); if (b<0.0)> { shh=-shh; } shh=c/(b+shh); } fg[0]=(sk+sm)*(sk-sm)-shh; fg[1]=sk*ek; for (i=kk; i<=mm-1; i++)=""> { sss(fg,cs); if (i!=kk) { e[i-2]=fg[0]; } fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1]; e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1]; fg[1]=cs[1]*s[i]; s[i]=cs[0]*s[i]; if ((cs[0]!=1.0)(cs[1]!=0.0)) { for (j=1; j<=n; j++)=""> { ix=(j-1)*n+i-1; iy=(j-1)*n+i; d=cs[0]*v[ix]+cs[1]*v[iy]; v[iy]=-cs[1]*v[ix]+cs[0]*v[iy]; v[ix]=d; } } sss(fg,cs); s[i-1]=fg[0]; fg[0]=cs[0]*e[i-1]+cs[1]*s[i]; s[i]=-cs[1]*e[i-1]+cs[0]*s[i]; fg[1]=cs[1]*e[i]; e[i]=cs[0]*e[i]; if (i { if ((cs[0]!=1.0)(cs[1]!=0.0)) { for (j=1; j<=m; j++)=""> { ix=(j-1)*m+i-1; iy=(j-1)*m+i; d=cs[0]*u[ix]+cs[1]*u[iy]; u[iy]=-cs[1]*u[ix]+cs[0]*u[iy]; u[ix]=d; } } } } e[mm-2]=fg[0]; it=it-1; } else { if (ks==mm) { kk=kk+1; fg[1]=e[mm-2]; e[mm-2]=0.0; for (ll=kk; ll<=mm-1; ll++)=""> { i=mm+kk-ll-1; fg[0]=s[i-1]; sss(fg,cs); s[i-1]=fg[0]; if (i!=kk) { fg[1]=-cs[1]*e[i-2]; &nbs v[ix]=d; } } } } else { kk=ks+1; fg[1]=e[kk-2]; e[kk-2]=0.0; for (i=kk; i<=mm; i++)=""> { fg[0]=s[i-1]; sss(fg,cs); s[i-1]=fg[0]; fg[1]=-cs[1]*e[i-1]; e[i-1]=cs[0]*e[i-1]; if ((cs[0]!=1.0)(cs[1]!=0.0)) { for (j=1; j<=m; j++)=""> { ix=(j-1)*m+i-1; iy=(j-1)*m+kk-2; d=cs[0]*u[ix]+cs[1]*u[iy]; u[iy]=-cs[1]*u[ix]+cs[0]*u[iy]; u[ix]=d; } } } } } } } return(1); } static void ppp(double a[],double e[],double s[],double v[],int m,int n) { int i,j,p,q; double d; if (m>=n) { i=n; } else { i=m; } for (j=1; j<=i-1; j++)="" {=""> a[(j-1)*n+j-1]=s[j-1]; a[(j-1)*n+j]=e[j-1]; } a[(i-1)*n+i-1]=s[i-1]; if (m { a[(i-1)*n+i]=e[i-1]; } for (i=1; i<=n-1; i++)="" {=""> for (j=i+1; j<=n; j++)=""> { p=(i-1)*n+j-1; q=(j-1)*n+i-1; d=v[p]; v[p]=v[q]; v[q]=d; } } return; } static void sss(double fg[2],double cs[2]) { double r,d; if ((fabs(fg[0])+fabs(fg[1]))==0.0) { cs[0]=1.0; cs[1]=0.0; d=0.0; } else { d=sqrt(fg[0]*fg[0]+fg[1]*fg[1]); if (fabs(fg[0])>fabs(fg[1])) { d=fabs(d); if (fg[0]<> { d=-d; } } if (fabs(fg[1])>=fabs(fg[0])) { d=fabs(d); if (fg[1]<0.0)> { d=-d; } } cs[0]=fg[0]/d; cs[1]=fg[1]/d; } r=1.0; if (fabs(fg[0])>fabs(fg[1])) { r=cs[1]; } else { if(cs[0]!=0.0) { r=1.0/cs[0]; } fg[0]=d; fg[1]=r; return; } } //复数矩阵相乘 //计算矩阵A(m*n)和B(n*k)的乘积，结果保存在C(m*k)中 //ar-长度为m*n的数组,存放A的实部 //ai-长度为m*n的数组,存放A的虚部 //br-长度为n*k的数组,存放B的实部 //bi-长度为n*k的数组,存放B的虚部 //cr-长度为m*k的数组，存放结果C的实部 //ci-长度为m*k的数组，存放结果C的虚部 void dbcmul(double ar[],double ai[],double br[],double bi[],int m,int n,int k,double cr[],double ci[]) { int i,j,l,u,v,w; double p,q,s; for (i=0; i<=m-1; i++)="" {=""> for (j=0; j<=k-1; j++)=""> { u=i*k+j; cr[u]=0.0; ci[u]=0.0; for (l=0; l<=n-1; l++)=""> { v=i*n+l; w=l*k+j; p=ar[v]*br[w]; q=ai[v]*bi[w]; s=(ar[v]+ai[v])*(br[w]+bi[w]); cr[u]=cr[u]+p-q; ci[u]=ci[u]+s-p-q; } } } return; } //复数矩阵求逆 //ar-长度为n*n的数组, n*n矩阵的实部 //ai-长度为n*n的数组, n*n矩阵的虚部 //n 矩阵的维数 int ddcinv(double ar[],double ai[],int n) { int *is,*js,i,j,k,l,u,v,w; double p,q,s,t,d,b; is=malloc(n*sizeof(int)); js=malloc(n*sizeof(int)); for (k=0; k<=n-1; k++)="" {=""> d=0.0; for (i=k; i<=n-1; i++)=""> { for (j=k; j<=n-1; j++)=""> { u=i*n+j; p=ar[u]*ar[u]+ai[u]*ai[u]; if (p>d) { d=p; is[k]=i; js[k]=j; } } } if (d+1.0==1.0) { free(is); free(js); printf("err**not inv\n"); return(0); } if (is[k]!=k) { for (j=0; j<=n-1; j++)=""> { u=k*n+j; v=is[k]*n+j; t=ar[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=ai[v]; ai[v]=t; } } if (js[k]!=k) { for (i=0; i<=n-1; i++)=""> { u=i*n+k; v=i*n+js[k]; t=ar[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=ai[v]; ai[v]=t; } } l=k*n+k; ar[l]=ar[l]/d; ai[l]=-ai[l]/d; for (j=0; j<=n-1; j++)=""> { if (j!=k) { u=k*n+j; p=ar[u]*ar[l]; q=ai[u]*ai[l]; s=(ar[u]+ai[u])*(ar[l]+ai[l]); ar[u]=p-q; ai[u]=s-p-q; } } for (i=0; i<=n-1; i++)=""> { if (i!=k) { v=i*n+k; for (j=0; j<=n-1; j++)=""> { if (j!=k) { u=k*n+j; w=i*n+j; p=ar[u]*ar[v]; q=ai[u]*ai[v]; s=(ar[u]+ai[u])*(ar[v]+ai[v]); t=p-q; b=s-p-q; ar[w]=ar[w]-t; ai[w]=ai[w]-b; } } } } for (i=0; i<=n-1; i++)=""> { if (i!=k) { u=i*n+k; p=ar[u]*ar[l]; q=ai[u]*ai[l]; s=(ar[u]+ai[u])*(ar[l]+ai[l]); ar[u]=q-p; ai[u]=p+q-s; } } } for (k=n-1; k>=0; k--) { if (js[k]!=k) { for (j=0; j<=n-1; j++)=""> { u=k*n+j; v=js[k]*n+j; t=ar[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=ai[v]; ai[v]=t; } } if (is[k]!=k) { for (i=0; i<=n-1; i++)=""> { u=i*n+k; v=i*n+is[k]; t=ar[u]; ar[u]=ar[v]; ar[v]=t; t=ai[u]; ai[u]=ai[v]; ai[v]=t; } } } free(is); free(js); return(1); } 四川大学计算机学院、软件学院 实 验 报 告 学生姓名 学号 班级 学生专业 计算机科学与技术 学院 课程名称 实 验 时 间 汇编语言程序设计(第七次) 实验项目 报告撰写时间 非数值运算程序设计 ? 掌握典型串操作指令的使用; 实验目的 ? 学习非数值运算程序设计。 ? X86系列桌面系统; 实验环境 ? UltraEdit-32、TASM、TLINK、TD。 ? 用stosw指令编写一字符显示模式下的清屏程序。 编程参考:DrctShw 实验内容 ? 用带重复前缀的串操作指令，完成教材219页的数据块传送程序，陈述源 数据区和目标数据区交叠情形如何避免数据丢失的方法;并对源程序作详 尽注释。 ? 字符模式清屏程序的源代码;、 stack0 segment stack 'stack' ; dw 40h dup(0) stack0 ends 实验记录 data segment Hello db "Hello,World!",0 Hi db "Hi! everybody.",0 1 db 1 x db 2 fctr equ 1218h y db 3 sum dw ? data ends code segment assume cs:code,ds:data,ss:stack0 main: mov ax,data mov ds,ax mov al,x add al,y mov sum,fctr mov dh,10 mov dl,8 mov ah,1eh lea si,Hello call DrctShw 2 mov ah,0 int 16h call ClrScr mov ah,0 int 16h mov dh,11 mov dl,8 mov ah,1eh lea si,Hi call DrctShw mov ah,0 int 16h mov ah,4ch int 21h DrctShw proc push bx push di push es 3 mov bx,0b800h mov es,bx mov bh,ah mov al,dh dec al mov bl,160 mul bl mov di,ax dec al mov al,dl mov bl,2 mul bl add di,ax mov ah,bh Next: mov al,[si] cmp al,0 jz rtrn mov es:[di],al 4 mov es:[di],ax inc di mov es:[di],bh inc di nc si jmp Next rtrn: pop es pop di . pop bx ret DrctShw endp ClrScr proc push ax push cx push di push es mov bx,0b800h mov es,bx 5 mov cx,80*25 mov di,0 mov ax,0020H cld Next0: rep stosw pop es pop di pop cx pop ax ret ClrScr endp code ends end main ? 数据传送程序的源代码; stck segment stack 'stack' dw 20 dup (0) stck ends 6 data segment srcblck db '128976598WS JKOPB:?65#RSTwVWI234' lgth equ $ - srcblck ;表示字符串srcblck的 blanck db lgth dup('$') ;操作后的字符串 crlf db 0dh,0ah,'$' ;执行换行的0dh和0ah的功能缓存区 data ends code segment assume cs:code,ds:data,ss:stck main: mov ax,data mov ds,ax mov es,ax lea dx,srcblck ;用9号功能显示源字符串 mov ah,9 int 21h lea dx,crlf ;用9号功能调用crlf字符串实现换行 mov ah,9 int 21h lea si,srcblck ;si指向源字符串首地址 7 lea di,srcblck+10 ;di指向操作字符串首地址 cld ;方向标志位清零 cmp si,di ;比较源字符串和操作字符串首地址大小 jae next ;若源地址大于操作地址则由首地址开始传送 Std ;否则，方向标志位置1，从串尾开始传送 add si,lgth-2 ;si指向源字符串串尾 add di,lgth-2 ;di指向操作字符串串尾 next:mov cx,lgth shr 1 ;定义rep重复次数 rep movsw ;重复传送串内容 lea dx,srcblck ;显示操作后串的内容 mov ah,9 int 21h mov ah,0 ;0号功能设置按键继续 int 16h mov ah,4ch int 21h code ends end main 8 ? 编程及调试中遇到的问题及解决办法。 清零时应考虑屏幕段基址以及屏幕共占的字节长度大小 出现传送串时方向失控，用cld将方向标志清零后恢复正常 反向传送串时串尾指向错误，将si.di指向si/di+lgth-2后正常 教师评议 成绩评定: 指导教师签名: 9 Matlab中巧用符号运算实现数值计算程序 的通用性 ISSN1008—9446 =『 第11卷第1期,2009年3月 Vo1.11,No.1,Mar.2009 Matlab中巧用符号运算实现数值计算 程序的通用性 刘兵,王辉,郝海燕 (承德石油高等专科学校数理系,河北承德,067000) 摘要:由于Matlab中内嵌了Maple的符号运算功能,因此在Matlab中实现数值计算程序的通用模式变得非常 方便.利用Matlab的符号运算功能,将求Jacobi矩阵(特别是涉及非线性函数)等繁杂的工作留给计算机完 成,从而大大简化了程序的规模,并实现了程序对不同函数的通用性. 关键词:Matlab;符号运算;数值计算 中图分类号:TP317文献标识码:B文章编号:1008—9446(2009)O1_0O55—03 UsingSymbolOperationtoMakeAvailableProgramsfor NumericalCalculationinMatlab LIUBing,WANGHui,HAOHai—yan (DepartmentofMathematics,ChengdePetroleumCollegeChengde067000,Hebei,China) Abstract:Forthetoolboxofsymbol,itispossibleforUStOmaketheavailableprogramsfor numericalcalculationinMatlab.Thispapersolvesthisproblembyusingthesymboloperation skillfully.Thuswenotonlypredigestthescaleoftheprogramsbutalsoachievethegoalofmaking programswhichispracticablefordifferentfunctions. Keywords:Matlab;symboloperation;numericalcalculation Matlab作为一门以矩阵为基本运算单位的数学软件,以其良好的交 互性和突出的矩阵运算能力成 为了数值计算领域中编程软件的首选,特别是其功能强大的工具箱 为用户提供了极大方便. 1预备知识 Matlab?中有许多的工具箱,每个工具箱中都有自己的专业性函数,符 号工具箱就是其中一种. Matlab中的符号工具箱源自于Maple,它拥有Maple中的大部分符号 运算功能.在数值计算中运用符号 运算,首先应掌握如下几点: 1)在Matlab中符号常量的定义形式为sym(x),其中X代表Matlab中 的符号常量. 2)符号变量的定义形式为y=sym(x),其中Y代表Matlab中的符号变量,该变量储存的内容为符 号常量,X. 3)函数feval:feval函数的表达形式为[a,b,…]=feval(f,x1,x2,…,xn),它的作用相当于我们在 Matlab的命令窗口中直接输入[a,b,…]=f(x1,x2,…,xn),表示求以(xl,x2,…,xn)为参数的函数f (x1,x2,…,xn)的返回值. 4)函数eval:eval函数表示形式为eval(f),其中f为字符串或符号表达式,它的作用为可以实现由 字符串或符号表达式到数值的转换.例如,在Matlab的命令窗口输入以下命令: 收稿日期:2008—10—20 作者简介:刘兵(1980一),男,河北隆化人,承德石油高等专科学校数理系,硕士研究生.研究方向:最优化理论及 其应用. ? 56-承德石油高等专科学校2009年第11卷第1期 >>Y=sym(X); >>f=y +5Y; >>X=3: >>eval(f) anS= 24 2Matlab中编写最优化通用程序 2.1一个求解无约束最优化问题的实例 下面,我们在Matlab平台上编写一个用牛顿法结合Armijo—Goldstein准则解一般无约束最优化问题 的通用程序模块.该程序只要求使用者输入将被优化的目标函数和初始点即可.操作过程如下: 1)借助符号运算功能求任何实值函数的梯度g以及hesse矩阵G的通用程序’GG.m’(模块1),具 体代码如下: function[g,G]=GG(f,X) n=length(X); fori=1:n d(i):sym([,Sint2str(i)],real,). end [W,e]=feval(f,X); l=jacobian(e,d); 1=1; G:jacobian(1,d); fori=1:n eva1([,Sint2str(i)=X(i)]); end g=eva1(1); G=eval(G); 2)用牛顿法结合Armijo.Goldstein准则解一般无约束最优化问题的主程序hiudun.m模块2),具 体代码如下: functionniudun(GG,f,X) [gG]=feval(GG,f,X); whilenorln(g)>=1e一10 a=1; b=10; P=1/4; while1 iffeval(f,x—ainv(G)g)<=feval(f,x)一Paginv(G)g ifferal(f,x—ainv(G)g)>=reval(f,x)一(1一P)$aginv(G)g t=a; display(t); break; e1se 刘兵,等:Matlab中巧用符号运算实现数值计算程序的通用性?57? a=(a+b)/2; continue; end else b=a; a=a/2: continue; end end ainv(d)g; X=X— [gG]=reval(GG,f,X); end disp(X); 说明:这里while循环部分为利用Armijo—Goldstein准则做线搜索的过程. 2.2数值实验 下面,通过计算一个具体函数来说明本文通用模块的应用.讨论目标函数 f(xI~72)=+:一1 ,初始点为(一1,0.5). 首先,在Matlab中编写关于函数f的一个In文件f.m,其中w代表f(x)的数值表达式,e代表f(x) 的符号表达式.为保证与前面模块的符号运算的一致性,此处函数返 回的函数符号表达式的常量仍为 s具体程序(模块3)代码为: function[W,e]=f(X) W=X(1)2+X(2)2一X(1)%X(2)/2; n=length(X); fori=1:n d(i)=sym([sint2str(i)],real); end e=d(1)-2+d(2)一d(1)d(2)/2; 注意,这里的符号变量数组d的每个分量d(i)对应值必为s. Matlab的命令窗口中输入niudun(‘GG’,’f’,[一1;0.5]),则看到如下结 果: >>niudun(GG,f,[一1;0.5]) S1= 一 l s2= 0.5000 t= 1 Sl= O s2= 0 (下转第60页) ? 60?承德石油高等专科学校2009年第11卷第1期 4)审核文件编制不规范.目前还没有对审核文件的要求,清洁生产审核机构往往缺少实地测量,为 了应付验收,把已经实施的工程列入清洁生产方案,而把真正需要而企业一时拿不出钱来实施的列入持 续清洁生产方案,完全是为审核而审核,为编制报告而编制报告,和企业一起欺骗管理部门,企业也就不 会认识到清洁生产审核的重要性,长此以往,造成管理部门责令企业进行清洁生产审核,企业同咨询服 务机构一起欺骗主管部门的恶性循环. 2对策及建议 2.1完善清洁生产审核法律制度 国家应在法律的角度强化清洁生产审核管理,明确主管部门,建立清洁生产审核服务机构和人员资 质认证体系,建立清洁生产审核相关的法律约束和处罚机制. 2.2加强清洁生产审核相关标准,导则等技术指导文件的制定与发布 当前,已有几十项清洁生产国家标准颁布,但仍不能满足日益发展的清洁生产审核工作的需要,应 加大标准的制定力度,同时应制定审核技术导则,定期发布重点行业的清洁生产技术规范,以使清洁生 产审核人员有据可依. 2.3加强宣传教育与培训 加强宣传教育与培训应该是全方位的教育,既包括对审核企业的教育,也包括对审核机构的培训, 同时还应包括对全社会,管理部门的教育培训,既包括清洁生产审核政策的宣传,又应该包括清洁生产 审核技术的培训,以此达到认识水平和技术水平同步提高. 2.4将清洁生产审核纳入到环境管理体系和能源管理体系中 结合现有环境管理和能源管理制度,研究清洁生产审核与这些制度相结合的全能性,为企业开展清 洁生产审核提供一定的优惠政策.将清洁生产审核纳入总量控制和节能减排工作之中,并逐步引入如 排污许可证制度等管理制度体系中去. 总之,清洁生产审核从方法学,技术科学等角度看都处于不断探索和完善中,是一项推动清洁生产 的重要工具和具体措施,解决好实施中存在的问题,对提升企业竞争力,保护环境,实现科学发展,可持 续发展必将起到更大的作用. 参考文献: [1]国家环境保护总局科技标准司.清洁生产审计培训教材[M].北京: 中国环境科学出版社,2001 [2]国家环境保护局.企业清洁生产审计手册[M].北京:中国环境科学出版社,1996. (上接第57页) O 0 显然,x=(0,0)就是我们所求的最优点. 由此看到,模块1和模块2就是通用模块,对任意二次可微函数,我们只需在模块3中把函数的具 体形式给出来即可,其他的象计算梯度和Hesse阵的工作就由通用模块完成,这是那些没有符号运算功 能的编程软件难以实现的.本文的技巧在于赋值语句d(i)=sym([‘s’int2str(i)],’real’),这是作者 花很长时间才摸索出来的,希望能为广大应用者提供方便. 参考文献: [1]葛哲学.精通MATLAB[M].北京.电子工业出版社.2008. [2]袁亚湘,孙文瑜.最优化理论与方法[M].北京.科学出版社.2005范文二：卷积的数值运算
范文三：[IT/计算机]数值计算程序大放送-矩阵运算
范文四：汇编语言第7次实验报告：非数值运算程序设计
范文五：[doc格式] Matlab中巧用符号运算实现数值计算程序的通用性
