b[j][i]=a[i][j];
}
return0.0;
}
ii.矩阵求逆的示例函数（C语言）
intinvGJ（double**a,intn）
{
int*is,*js,i,j,k,l,u,v;
doubled,p;
is=（int*）malloc（n*sizeof（int））;
js=（int*）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[i][j]）;
if（p>d）
{d=p;is[k]=i;js[k]=j;}
}
if（d+1.0==1.0）
{
free（is）;free（js）;printf（"errornotinv\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[k][j];a[k][j]=a[is[k]][j];a[is[k]][j]=p;
}
if（js[k]!
=k）
for（i=0;i<=n-1;i++）
{
u=i*n+k;v=i*n+js[k];
p=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=p;
}
l=k*n+k;
a[k][k]=1.0/a[k][k];
for（j=0;j<=n-1;j++）
if（j!
=k）
{
u=k*n+j;a[k][j]=a[k][j]*a[k][k];
}
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[i][j]=a[i][j]-a[i][k]*a[k][j];
}
for（i=0;i<=n-1;i++）
if（i!
=k）
{
u=i*n+k;a[i][k]=-a[i][k]*a[k][k];
}
}
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[k][j];a[k][j]=a[js[k]][j];a[js[k]][j]=p;
}
if（is[k]!
=k）
for（i=0;i<=n-1;i++）
{
u=i*n+k;v=i*n+is[k];
p=a[i][k];a[i][k]=a[i][is[k]];a[i][is[k]]=p;
}
}
free（is）;free（js）;
return
（1）;
}intinvGJ（double**a,intn）
{
int*is,*js,i,j,k,l,u,v;
doubled,p;
is=（int*）malloc（n*sizeof（int））;
js=（int*）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[i][j]）;
if（p>d）
{d=p;is[k]=i;js[k]=j;}
}
if（d+1.0==1.0）
{
free（is）;free（js）;printf（"errornotinv\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[k][j];a[k][j]=a[is[k]][j];a[is[k]][j]=p;
}
if（js[k]!
=k）
for（i=0;i<=n-1;i++）
{
u=i*n+k;v=i*n+js[k];
p=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=p;
}
l=k*n+k;
a[k][k]=1.0/a[k][k];
for（j=0;j<=n-1;j++）
if（j!
=k）
{
u=k*n+j;a[k][j]=a[k][j]*a[k][k];
}
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[i][j]=a[i][j]-a[i][k]*a[k][j];
}
for（i=0;i<=n-1;i++）
if（i!
=k）
{
u=i*n+k;a[i][k]=-a[i][k]*a[k][k];
}
}
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[k][j];a[k][j]=a[js[k]][j];a[js[k]][j]=p;
}
if（is[k]!
=k）
for（i=0;i<=n-1;i++）
{
u=i*n+k;v=i*n+is[k];
p=a[i][k];a[i][k]=a[i][is[k]];a[i][is[k]]=p;
}
}
free（is）;free（js）;
return
（1）;
}
iii.矩阵求逆函数的调用（C语言）
#include
#include
#include
intinvGJ（double**a,intn）;
voidmain（）
{
inti,j;
double**AA;
//首先对二维指针Naa分配内存，采用C语言的方法
/*AA=（double**）malloc（sizeof（double）*2）;
for（i=0;i<2;i++）
{
AA[i]=（double*）mallo（sizeof（double）*2）;
}
*/
//首先对二维指针Naa分配内存，采用C++语言的方法
AA=newdouble*[2];
for（i=0;i<2;i++）
{
AA[i]=newdouble[2];
}
doubleBB[2][2]={1,2,3,4};
for（i=0;i<2;i++）
{
for（j=0;j<2;j++）
{
AA[i][j]=BB[i][j];
}
}
//调用矩阵求逆函数
invGJ（AA,2）;
printf（"矩阵AA的逆阵如下\n"）;
for（i=0;i<2;i++）
{
for（j=0;j<2;j++）
{
printf（"%10.4lf",AA[i][j]）;
}
printf（"\n"）;
}
doubleCC[2][2];
printf（"AA与其逆阵的乘积如下（理论上是单位阵）\n"）;
for（i=0;i<2;i++）
{
for（j=0;j<2;j++）
{
CC[i][j]=0.0;
for（intk=0;k<2;k++）
{
CC[i][j]+=AA[i][k]*BB[k][j];
}
printf（"%10.4lf",CC[i][j]）;
}
printf（"\n"）;
}
//C语言释放AA二维指针的方法
/*for（i=0;i<2;i++）
{
free（AA[i]）;
}
free（AA）;
*/
//C++语言释放AA二维指针的方法
for（i=0;i<2;i++）
{
deleteAA[i];
}
deleteAA;
}
实验3误差椭圆元素计算
一、实验名称：
误差椭圆元素计算。
二、实验目的和任务：
掌握误差椭圆和相对误差椭圆元素的计算公式，并采用C或者C++语言变成实现。
三、实验要求：
1每人独立编写出误差椭圆和相对误差椭圆元素的计算程序，并调试通过；
2采用VC++6.0开发平台，C或者C++语言编写程序；
3写出计算的结果。
4本实验属于综合性，设计性实验，对学生的要求比较高，会综合使用矩阵加法，乘法以及转置和求逆的程序。
5注意事项：
坐标方位角计算时必须考虑X与Y所在的象限。
四、实验内容：
参考教材《测量平差》P103-P112。
具体的例子可以参考P111页的例题，进行程序的调试。
（一）误差椭圆元素计算公式：
（二）相对误差椭圆元素计算公式：
五、例子（P111）
在某三角网中插入P1及P2两个新点。
设用间接平差法平差。
平差之后这两点之间的协因数阵如下：
根据以上的公式，分别计算出未知点P1和P2的误差椭圆元素以及这两点之间的相对误差椭圆元素。
实验4水准网间接平差程序设计
一、实验名称：
水准网间接平差程序设计
二、实验目的：
掌握用间接平差法对任意网形的水准网进行平差的算法设计以及程序编制；并学习采用读文件处理数据的方法。
三、实验任务：
用C/C++编写水准网间接平差程序，并调试通过，用实测数据实算分析；同时评定精度。
四、实验要求：
1个人独立编写程序，原始数据存放于文本文件或者数据库文件中；
2程序应该具有通用性，即任意网形都可以平差；
3用实测的数据进行计算分析。
4评定待定点的高程精度。
五、实验内容
1数据文件的编制格式
总点数未知点数测段数（观测高差的个数）
已知点1点名高程
已知点2点名高程
未知点3点名0
未知点4点名0
未知点5点名0
。
。
。
起点点号终点点号观测高差路线长度（km）
起点点号终点点号观测高差路线长度（km）
起点点号终点点号观测高差路线长度（km）
。
。
。
2平差原理
由观测值的起始和终点号或者水准网网形，形成误差方程的系数矩阵B（也叫设计矩阵），由观测的路线长度形成观测高差的权阵P（观测值独立，P为对角阵），
原理如下：
（C为任意常数）
3法方程系数矩阵与闭合差的自动累加
由于N和W具有可加性，即每读取一个高差观测值，即可得到一个误差方程的系数向量，然后累加到法方程系数矩阵与闭合差中；当高差观测值读取完毕的时候，法方程系数矩阵与闭合差也累加完毕；此时可用公式直接计算待定点高程的改正数，加上高程近似值，就得到了高程的平差值。
4精度评定
验后单位权中误差：
待定点高程（未知参数）的协因数阵：
待定点高程的方差：
改正数V的协因数阵：
高差平差值的协因数阵：
5例子1（P79）
（1）in.txt文件
435
1A237.483
2B0.0
3C0.0
4D0.0
125.8353.5
233.7822.7
139.6404.0
437.3843.0
142.2702.5
（2）源程序
//Gckzwpc.h头文件代码
classCGckzwpc
{
public:
CGckzwpc（）;
virtual~CGckzwpc（）;
public:
boolReadData（CStringfilename）;//读水准网平差数据文件
voidpc（）;//平差函数
voidjdpd（）;//精度评定函数
intinvGJ（double**a,intn）;//求逆函数
intnz,nw,ne,nn;//文件头信息：
总点数，未知点数，已知点数，测段数
intn1[20],n2[20];//存放高差起点与终点的点号
doubleH[50],h[50],W[20],X[20],B[50][20],V[50],VPV;
double**Nbb,S[50],l[50],P[50],ph[50];
CStringdm[20];//控制点点名
doubleSIG0,DX[20][20];
};
//Gckzwpc.cpp代码
boolCGckzwpc:
:
ReadData（CStringfilename）
{
inti;
intMAXLINE=512;
charbuff[513],ch1[15];
CStdioFilefp;
if（!
fp.Open（filename,CFile:
:
modeRead|CFile:
:
typeText,NULL））
{
AfxGetApp（）->m_pMainWnd->MessageBox（
"数据文件不存在或数据文件错!
",
"进程......!
!
!
",MB_OK|MB_ICONSTOP）;
_exit
（1）;
returnFALSE;
}
fp.ReadString（buff,MAXLINE）;
sscanf（buff,"%d%d%d",&nz,&nw,&nn）;
ne=nz-nw;
intdh;
doublegc;
for（i=0;i{
fp.ReadString（buff,MAXLINE）;
sscanf（buff,"%d%s%lf",&dh,ch1,&gc）;
dm[i]=ch1;H[i]=gc;
}
for（i=0;i{
fp.ReadString（buff,MAXLINE）;
sscanf（buff,"%d%d%lf%lf",
&n1[i],&n2[i],&h[i],&S[i]）;
}
fp.Close（）;
returnTRUE;
}
voidCGckzwpc:
:
pc（）
{
inti,j,k;
//doubleNbb1[20][20];
Nbb=newdouble*[nw];
for（i=0;iNbb[i]=newdouble[nw];
for（i=0;i{
P[i]=0.0;
l[i]=0.0;
for（j=0;jB[i][j]=0.0;
}
for（i=0;i{
if（n1[i]>ne）B[i][n1[i]-ne-1]=-1;
if（n2[i]>ne）B[i][n2[i]-ne-1]=+1;
P[i]=1/S[i];
l[i]=H[n1[i]-1]+h[i]-H[n2[i]-1];
}
for（i=0;i{
for（j=0;j{
Nbb[i][j]=0.0;
for（k=0;k{
Nbb[i][j]+=B[k][i]*B[k][j]*P[k];
}
}
}
for（i=0;i{
W[i]=0.0;
for（k=0;k{
W[i]+=B[k][i]*P[k]*l[k];
}
}
invGJ（Nbb,nw）;
for（i=0;i{
X[i]=0.0;
for（k=0;k{
X[i]+=Nbb[i][k]*W[k];
}
}
}
//精度评定函数
voidCGckzwpc:
:
jdpd（）
{
inti,j,k;
FILE*fp;
fp=fopen（"out.txt","w"）;
VPV=0.0;
for（i=0;i{
V[i]=0.0;
for（k=0;k{
V[i]+=B[i][k]*X[k];
}
V[i]+=-l[i];
VPV+=V[i]*V[i]*P[i];
ph[i]=h[i]+V[i];
}
//计算验后单位权中误差
SIG0=sqrt（VPV/（nn-nw））;
//计算未知参数的方差
for（i=0;i{
for（j=0;j{
DX[i][j]=SIG0*SIG0*Nbb[i][j];
}
}
doubleQLL[50][50],SIGL[50][50];
//计算观测值的平差值的中误差QLL=B*inv（Nbb）*BTDLL=SIG0*SIG0*QLL
for（i=0;i{
for（j=0;j{
QLL[i][j]=0.0;
SIGL[i][j]=0.0;
for（k=0;k{
for（intm=0;m{
QLL[i][j]+=Nbb[k][m]*B[i][k]*B[j][m];
}
}
SIGL[i][j]=SIG0*sqrt（QLL[i][j]）;
}
}
fprintf（fp,"********水准网间接平差结果**********\n\n"）;
fprintf（fp,"总点数＝%3d未知点数＝%3d测段数＝%3d\n",
nz,nw,nn）;
fprintf（fp,"\n验后单位权中误差＝%6.2lf（mm）\n\n",SIG0*1000）;
fprintf（fp,"\n起点号终点号观测高差（m）路线长（km）改正数（mm）平差高差（m）中误差（mm）\n"）;
for（i=0;i{
fprintf（fp,"%3d%3d%8.4lf%8.4lf%8.2lf%8.4lf%8.2lf\n",
n1[i],n2[i],h[i],S[i],V[i]*1000,ph[i],SIGL[i][i]*1000）;
}
fprintf（fp,"\n\n已知点已知高程（m）\n"）;
for（i=0;i
升级会员 