return (i >= j) ? i*(i + 1) / 2 + j : j*(j + 1) / 2 + i; } ////////////////////////////////////////////////////// //显示提示信息 void MyBreak(char *fmt, ...) {
Q[i] = P[i];
if (inverse(Q, n) == false)
5
Jack 专属
{ MyBreak("观测值的权矩阵奇异！"); return -1.0;
}
for (int i = 0; i < r; i++) {
for (int j = 0; j <= i; j++) {
double nij = 0.0; for (int k = 0; k < n; k++) {
}
fprintf(fp, fmt, A[i]); } fprintf(fp, "\n"); } ///////////////////////////////////////////////////////////////////// // 向文件输出对称矩阵（下三角元素） void PrintM2(FILE *fp, double M[], int n, int t, char *fmt, char *title, bool IsLabel) { if (title)
//法方程系数矩阵求逆 if (inverse(N, t) == false) {
delete[]U; return -1.0; }
//计算X for (int i = 0; i < t; i++)
4
Jack 专属
{ X[i] = 0.0; for (int j = 0; j < t; j++) X[i] += N[ij(i, j)] * U[j];
#endif//VC_EXTRALEN
}
/////////////////////////////////////////////// // 对称正定矩阵求逆（仅存下三角元素） bool inverse(double a[], int n) {
double *a0 = new double[n]; for (int k = 0; k < n; k++) {
fprintf(fp, "\n%s:", title);
int index = 0; for (int i = 0; i < n; i++) {
if (IsLabel) fprintf(fp, "\n%3d", i + 1);
else fprintf(fp, "\n ");
for (int j = 0; j <= i; j++) {
if (inverse(N, r) == false) {
delete[]Q; delete[]N;
6
return -1.0; }
Jack 专属
double m = Co_Adjust(r, t, A, W, N, X, QX, V); for (int i = 0; i < t; i++)
X[i] = -X[i];
char buffer[256]; va_list argptr; va_start(argptr, fmt); vsprintf_s(buffer, fmt, argptr);def VC_EXTRALEN AfxMessageBox(buffer);
#else printf(buffer); getchar();
double pvv = Calculate_q(P, V, n);//权倒数计算
return sqrt(pvv / (n - t));
}
///////////////////////////////////////////////////////////// // 具有参数的条件平差 double Condition_With_X(int n, int r, int t, bool correlative, double B[], double A[], double W[], double P[], double X[], double QX[], double V[]) {
ki += N[ij(i, j)] * V[j]; } K[i] = ki; }
//计算向量V double pvv = 0.0; if (correlative)//观测值相关 {
for (int i = 0; i < n; i++) {
double vi = 0.0; for (int k = 0; k < n; k++) {
if (j>0 && j%t == 0) fprintf(fp, "\n ");
fprintf(fp, fmt, M[index++]); }
} fprintf(fp, "\n"); }
////////////////////////////////////////////////////////// // 权倒数计算函数 double Calculate_q(double *Q, double *F, int t) {
for (int s = 0; s < r; s++) {
vi += Q[ij(i, k)] * B[s*n + k] * K[s];
7
Jack 专属
} } V[i] = vi; } pvv = Calculate_q(P, V, n);
} else //观测值独立 {
for (int i = 0; i < n; i++) {
else a0[i] = ai0 / a00;
for (int j = 1; j <= i; j++) {
a[(i - 1)*i / 2 + j - 1] = a[i*(i + 1) / 2 + j] + ai0*a0[j];
} } for (int i = 1; i < n; i++) {
a[(n - 1)*n / 2 + i - 1] = a0[i]; } a[n*(n + 1) / 2 - 1] = 1.0 / a00; } delete[]a0; return true; }
/////////////////////////////////////////////////////////////////////////////// //具有参数的条件平差算例 void Condition_With_X_Example() {
FILE *fp; fopen_s(&fp, "data.txt", "r"); if (fp == NULL) {
double *Q = new double[n*(n + 1) / 2];//观测值的权逆阵 double *N = new double[r*(r + 1) / 2];//法方程系数矩阵
//
计算矩阵N
if (correlative)//观测值相关
{
for (int i = 0; i < n*(n + 1) / 2; i++)
for (int s = 0; s < n; s++) {
for (int i = 0; i < t; i++) {
U[i] += A[k*t + i] * P[ij(k, s)] * L[s]; for (int j = 0; j <= i; j++) {
N[ij(i, j)] += A[k*t + i] * P[ij(k, s)] * A[s*t + j]; } } } }
if (title) fprintf(fp, "\n %s: ", title);
int j = 0; for (int i = 0; i < size; i++) {
if (i%t == 0) {
j++; if (IsLabel)
fprintf(fp, "\n%3d", j);
2
Jack 专属
else fprintf(fp, "\n");
double a00 = a[0]; if (a00 + 1.0 == 1.0)
1
Jack 专属
{ delete[]a0; return false;
} for (int i = 1; i < n; i++) {
double ai0 = a[i*(i + 1) / 2];
if (i <= n - k - 1) a0[i] = -ai0 / a00;
if (m < 0.0) {
delete[]Q; delete[]N; return -1.0;
}
double *K = new double[r]; for (int i = 0; i < r; i++) {
double ki = 0.0; for (int j = 0; j < r; j++) {