Алгоритмическая и программная реализация прогноза почвенных параметров

При числе параметров более трех и порядка полиномов а_S(X_S) более второго уравнение (3.1.2) становится весьма громоздким и малоэффективным в использовании, оно излишне перегружено суммой произведений различных комбинаций факторов. Поэтому был разработан упрощенный метод описания многофакторных процессов.

3.2 Алгоритм метода

Точность математического описания в конечном итоге определяет коэффициент детерминации. Приемлемым может быть тот алгоритм описания, который позволяет получить аппроксимирующую функцию процесса с удовлетворительной предсказательной способностью (т. е. если коэффициент детерминации более 0,8).

Более удобным в данном случае может оказаться метод так называемого “последовательного разглаживания гиперповерхности отклика”. Идея этого метода заключается в последовательном снятии эффектов всех переменных полиномиального уравнения (3.1.1).

В основу положен ступенчатый регрессионный метод [3] при множественной линейной регрессии.

В этом случае уравнение (3.1.1) можно записать в виде:

Y=сt_j(X_j)(3.2.1)

где t_j(X_j) -- полиномы возрастающих степеней переменной X_j.

Перед вычислением полиномов t_j необходимо проводить предварительное ранжирование переменных по степени влияния их на параметр оптимизации y. Начинать вычисление полинома нужно для наиболее сильно влияющей переменной. Расчет коэффициентов полинома и оптимальный порядок его целесообразно производить, используя ортогональные многочлены Чебышева .

Путь задан процесс в виде матрицы наблюдений:

y₁ X₁₁ . . . X_j1 . . . X_m1

. . . . . . . . . . . . . .

y_n X_1n . . . X_jn . . . X_mn

где m -- число независимых параметров;

n -- число наблюдений;

y -- зависимый параметр;

Х -- матрица наблюдений независимых параметров.

1. Последовательно для всех пар массивов y-X₁, y-X₂, . . .y-X_m вычислить остаточную дисперсию D_(k-1)1, D_(k-1)2, . . ., D_(k-1)m, характеризующую степень влияния независимого параметра на зависимый y.

Для этого воcпользуемся ортогональными многочленами Чебышева.

Искомая зависимость условного среднего y от Х_j имеет вид:

y_k=t_j(X_j)=сb_a*S_a(X_j) (3.2.2)

где S_k(X_j) -- ортогональный многочлен Чебышева k-го порядка.

S₀(X_j)=1 -- полином нулевой степени.

Общая дисперсия:

D₀=1/nс(y_i-b₀)²(3.2.3)

где

b₀=1/nсy_i, S₁(X_j)=X_j+f₁,

f₁=-1/nсX_ji

b₁=сy_i*S₁(X_ji)/с[S₁(X_ji)]², y_k=b₀*S₀(X_j)+b₁*S₁(X_j)

Остаточная дисперсия:

D_kj=1/nс(y_i-y_ki)²,

где k=1.

Если D_(k-1)j/D_kj<=F_kp, где F_kp -- табличное значение для распределения Фишера, то осуществляется переход к следующему j.

Иначе k=k+1 и переход к пункту 2.

2. Вычисляем:

S_k(X_j)=(X_j+B_k)*S_k-1(X_j)-g_k*S_k-2(X_j),

где

B_k=-сX_ji[S_k-1(X_ji)]²/с[S_k-1(X_ji)]², (3.2.4)

g_k=сX_ji*S_k-2(X_ji)*S_k-1(X_ji)/с[S_k-2(X_ji)]²,

y_k=y_k-1+b_k*S_k(X_j),

b_k=сy_i*S_k(X_ji)/с[S_k(X_ji)]²,

Если D_(k-1)j/D_kj<=F_kp, то запоминается D_(k-1)j, иначе k=k+1 и возврат к пункту 2 и так до тех пор, пока D_(k-1)j<=F_kpD_kj для всех j=1,2,...,m.

3. Выбрать из всех D_(k-1)j минимальное, где j=1,2,...,m.

Если минимальных остаточных дисперсий получится несколько равных между собой, то выбрать последнюю при переборе.

4. Если D_(k-1)j=D_n,начальной промежуточной дисперсии, то конец решения. Иначе перейти к пункту 5.

5. Произвести как бы разглаживание поверхностей отклика в направлении переменной Х_j, вычитая из выборочных значений y_i величины, рассчитанные по t_j(X_j).

Сформировать массив:

y_1i=y_i-t_j(X_j)

6. Заменить в матрице наблюдений массив y_i на y_1i.

7. Повторить вычисления с пункта 1 до пункта 6 для матрицы наблюдений с учетом замененных массивов y_i на y_1i до тех пор, пока последняя остаточная дисперсия D_(k-1)j=D/r уменьшится в r раз по сравнению с начальной общей дисперсией величины y.

Такая процедура продолжается до тех пор, пока не будет осуществлено элиминирование всех переменных или получено уравнение с заданной предсказательной способностью.

3.3 Определение вкладов

Согласно полиномиальному приближению (3.2.3) имеем для зависимой переменной:

y=с сb_aj*S_a(X_j)(3.3.1)

где a -- степень полинома;

j -- номер независимой переменной;

K_j -- максимальная степень полинома для j--ой переменной;

b_aj -- коэффициент разложения.

Для каждого наблюдения по (3.3.1) находим значение зависимой переменной:

y_i=с сb_aj*S_a(X_ji), i=1,2, ... ,n.(3.3.2)

где n -- число наблюдений.

Определим сумму всех значений зависимой переменной:

yS=сy_i(3.3.3)

Аналогичные суммы находим для независимых переменных:

XS_j=с сb_aj*S_a(X_ji), где j=1,2,...,m.(3.3.4)

Из (3.3.3) и (3.3.4) получаем вектор вкладов для независимых переменных в полиномиальной модели (3.3.1):

W_j=XS_j/YS(3.3.5)

Вектор W дает физическую оценку вклада каждой независимой переменной в полиномиальной модели для y.

Он дает возможность определить степень обусловленности зависимой переменной от независимой.

3.4 Описание исходных данных

АХ(J,I)-общая матрица наблюдений, где J-номер параметра, I-номер наблюдения,

Х(J,I)-матрица независимых параметров,

Y(I)-вектор значений зависимого параметра,

MNN(I)-вектор номеров независимых параметров,

INOY-номер зависимого параметра,

ZZZ-вспомогательный параметр, равный 1.

ККК-максимальная степень аппроксимирующих полиномов,

КS1-максимальная степень аппроксимирующих полиномов для расчета вкладов (обычно ККК = КS1)

FKP-значение параметра в критерии Фишера (обычно 1),

FSS-верхняя граница для отношения D₀/D_кон.

Для хорошей аппроксимации FSS=10.

N-число реализаций.

3.5 Описание переменных в программе(приложение-1)

АК0(I)-массив b₀ для каждой независимой переменной, I=1,2,...,М.

АКК(I,J)-массив b_ij для i-ой независимой переменной и для j-ой степени полинома; i=1,2,...,М. j=1,2,...20.

А00-начальное среднее зависимой переменной Y(I)

DN0-начальная дисперсия зависимой переменной Y(I)

N-число реализаций,

М-число независимых переменных,

А0-текущее начальное среднее зависимой переменной Y(I),b₀

DN-текущая начальная дисперсия зависимой переменной

ALF1-f₁.

S(1,I)-I-ое значение полинома первой степени,

АСН-числитель дроби, b₁.

AZ-знаменатель дроби для b₁.

АI-коэффициент b₁.

К-порядок полинома.

YK(I)-массив модельных значений зависимой переменной, полученный при аппроксимации по одной независимой переменной I=1,2,...,N.

D(J)-текущее значение дисперсии между Y(I) и YK(I), J=1,2...,M.

DТ-предыдущее значение текущей дисперсии D(J).

ВКС-числитель дроби для Вк.

ВКZ-знаменатель дроби для Вк.

ВК-Вк.

GC-числитель для дроби g_k.

GZ-знаменатель дроби для g_к.

GK-g_k.

АК-числитель дроби для b_k, k=2,3,...,

АКZ-знаменатель дроби для b_k, k=2,3,...,

АК-b_k ,k=2,3,...,

ILS-программный переключатель.

XMI-минимальная дисперсия из D(J), J=1,2,...,М.

KS-номер минимальной дисперсии

АКМ(1)-b₁

ALM1(J)-массив f₁,J=1,2,...,М.

ВС-Вк

GC-g_k

ВСМ1(J,КК)-значения В_j,kk для j-ой независимой переменной для полинома КК-ого порядка

GCM1(J,KK)-значения g_j,kk для j-ой независимой переменной для полинома КК-ого порядка

ВВ(I)-рабочий вектор, I=1,2,...,N.

В(I)-рабочий вектор, I=1,2,...,N.

SM1(J)-рабочий вектор, J=1,2,...,M.

SM2(J)-вектор вкладов, J=1,2,...,M.

YZK(I)-рабочий вектор, I=1,2,...,N.

3.6 Программная реализация упрощенного метода Брандона

Для реализации метода Брандона была написана подпрограмма на алгоритмическом языке Фортран--V и внешняя вызываемая функция на языке С++.

Текст подпрограммы представлен в приложении-1.

Обращение: CALL MNLE (ZZZ,KKK,KS1,FSS,N,KJ,FKP,AX,MNN,INOY).

Описание параметров (входных):

ZZZ-параметр, равный 1.

ККК-задаваемая максимальная степень аппроксимирующих полиномов,

КS1-задаваемая максимальная степень аппроксимирующих полиномов для расчета вкладов (обычно ККК = КS1)

FSS-верхняя граница для отношения D₀/D_кон. Для хорошей аппроксимации FSS>10.

N-число реализаций,

KJ-число независимых переменных,

FKP-значение параметра в критерии Фишера

АХ(J,I)-общая матрица наблюдений,

где J-номер параметра, I-номер реализации

MNN(I)-вектор номеров независимых параметров размера КJ.

INOY-номер зависимого параметра.

Используемые подпрограммы: нет.

1.Проанализированы методы построения полиномиальных моделей аппроксимации.

2.Показана невыгодность использования полиномиальных моделей с большим количеством параметров-аргументов из-за больших вычислительных ошибок,связанных с конечностью разрядной сетки ЭВМ.

3.Описан алгоритм упрощенного метода Брандона Д. с автоматическим выбором степеней аппроксимирующих полиномов.

4.Приводится разработанный в исследовании алгоритм для определения вкладов параметров-аргументов в полиномиальных моделях.

5.Дано описание входных и выходных параметров подпрограммы по реализации упрощенного метода Брандона Д.

Глава 4. Результаты анализа

4.1 Анализ парных и групповых обусловленностей при временном изменении почвенных параметров

Для анализа были взяты значения почвенных параметров по пахотному слою(глубина взятия образца 0-27см) для чернозема южного малогумусного среднемощного тяжелосуглинистого в 1983 году и значения почвенных параметров по пахотному слою(глубина взятия образца 0-22см) для чернозема южного ср/мощного глинистого в 1963 году совхоза Свердлова Тоцкого района Оренбургской области.

Значения рассматриваемых почвенных параметров приведены в таблице-1.

Таблица-1

Параметр	Значение в 1963г	Значение в 1983г
Гумус по Тюрину	6,84	4,6
Поглощенный Са	20,74	18
Поглощенный Мg	7,88	2,4
Содержание фракций 1-0,25	0,72	1,1
Содержание фракций 0,25-0,05	9,63	11,4
Содержание фракций 0,05-0,01	25,72	30,7
Содержание фракций 0,01-0,005	10,89	4,3
Содержание фракций 0,005-0,001	20,12	18,3
Содержание фракций менее 0,001	39,92	34,2

Были построены две нормализованные матрицы исследования со следующими характеристиками нормальных распределений[1]:

Таблица-2

Параметр	Значение математического ожидания в 1963 году	Значение среднего квадратического отклонения в 1963 году	Значение математического ожидания в 1983 году	Значение среднего квадратического Отклонения в 1983 году
Гумус по Тюрину	6,84	0,1	4,6	0,01
Поглощенный Са	20,74	1,1	18	1,2
Поглощенный Мg	7,88	1,1	2,4	0,7
Содержание фракций 1-0,25	0,72	0,02	1,1	0,02
Содержание фракций 0,25-0,05	9,63	0,12	11,4	0,3
Содержание фракций 0,05-0,01	25,72	0,1	30,7	0,5
Содержание фракций 0,01-0,005	10,89	0,1	4,3	0,02
Содержание фракций 0,005-0,001	20,12	0,1	18,3	0,02
Содержание фракций менее 0,001	39,92	1,5	34,2	0,6

Для определения парных обусловленностей на построенных матрицах исследования был проведен корреляционный анализ.

Результаты корреляционного анализа на матрице 1963 года:

параметр 1-(Гумус по Тюрину в %)

с параметром-(Гумус по Тюрину в %) коэффициент корреляции= 1.000

Результаты корреляционного анализа на матрице 1983 года:

параметр 1-(Гумус по Тюрину в %)

с параметром-(Гумус по Тюрину в %) коэффициент корреляции= 1.000

Как видно из результатов корреляционного анализа ,гумус не имеет даже слабой парной обусловленности, при которой коэффициент парной корреляции r лежит на отрезке [0,3;0,5].

Для определения групповых обусловленностей на нормализованных матрицах исследования был проведен факторный анализ.

Результаты факторного анализа на матрице 1963 года:

Таблица-3-Объединение по фактору 5

в факторе- 5 базовый параметр- 1 (Гумус по Тюрину в %)

Таблица-4-Объединение по фактору 3

в факторе- 3 базовый параметр- 2 (Поглощенный кальций в мг-экв)

Таблица-5-Объединение по фактору 6

в факторе- 6 базовый параметр- 3 (Поглощенный магний в мг-экв)

Таблица-6-Объединение по фактору 1

в факторе- 1 базовый параметр- 7 (Фракция 0,01-0,005 мм в %)

Таблица-7-Объединение по фактору 4

в факторе- 4 базовый параметр- 5 (Фракция 0,25-0,05 мм в %)

Таблица-8-Объединение по фактору 7

в факторе- 7 базовый параметр- 9 (Фракция менее 0,001 мм в %)

Таблица-9-Объединение по фактору 2

в факторе- 2 базовый параметр- 8 (Фракция 0,005-0,001 мм в %)

Результаты факторного анализа на матрице 1983 года:

Таблица-10-Объединение по фактору 5

в факторе- 5 базовый параметр- 1 (Гумус по Тюрину в %)

Таблица-11-Объединение по фактору 2

в факторе- 2 базовый параметр- 2 (Поглощенный кальций в мг-экв)

Таблица-12-Объединение по фактору 6

в факторе- 6 базовый параметр- 3 (Поглощенный магний в мг-экв)

Таблица-13-Объединение по фактору 1

в факторе- 1 базовый параметр- 7 (Фракция 0,01-0,005 мм в %)

Таблица-14-Объединение по фактору 4

в факторе- 4 базовый параметр- 5 (Фракция 0,25-0,05 мм в %)

Таблица-15-Объединение по фактору 7

в факторе- 7 базовый параметр- 9 (Фракция менее 0,001 мм в %)

Таблица-16-Объединение по фактору 3

в факторе- 3 базовый параметр- 8 (Фракция 0,005-0,001 мм в %)

Факторный анализ показал отсутствие больших объединений по факторам и для матрицы 1963 года и для матрицы 1983 года.

То есть факторный анализ показал отсутствие групповых обусловленностей между почвенными параметрами.

Ранжируем по возрастанию элементы столбцов нормализованных матриц исследования и проведем на них корреляционный и факторный анализы.

Результаты предварительной обработки и корреляционного анализа на ранжированной матрице исследования по 1963 году:

Чернозем южный ср/мощный глинистый, глубина 0-22см, 1963 год

Результаты предварительной обработки



Корреляционная матрица R.

параметр 1-(Гумус по Тюрину в %)

с параметром-(Гумус по Тюрину в %)

коэффициент корреляции= 1.000

с параметром-(Поглощенный кальций в мг-экв)

коэффициент корреляции= .966

с параметром-(Поглощенный магний в мг-экв)

коэффициент корреляции= .960

с параметром-(Фракция 1-0,25 мм в %)

коэффициент корреляции= .980

с параметром-(Фракция 0,25-0,05 мм в %)

коэффициент корреляции= .957

с параметром-(Фракция 0,05-0,01 мм в %)

коэффициент корреляции= .985

с параметром-(Фракция 0,01-0,005 мм в %)

коэффициент корреляции= .970

с параметром-(Фракция 0,005-0,001 мм в %)

коэффициент корреляции= .978

с параметром-(Фракция менее 0,001 мм в %)

коэффициент корреляции= .983

Результаты корреляционного анализа на ранжированной матрице исследования по 1983 году:

параметр 1-(Гумус по Тюрину в %)

с параметром-(Гумус по Тюрину в %)

коэффициент корреляции= 1.000

с параметром-(Поглощенный кальций в мг-экв)

коэффициент корреляции= .965

с параметром-(Поглощенный магний в мг-экв)

коэффициент корреляции= .961

с параметром-(Фракция 1-0,25 мм в %)

коэффициент корреляции= .981

с параметром-(Фракция 0,25-0,05 мм в %)

коэффициент корреляции= .958

с параметром-(Фракция 0,05-0,01 мм в %)

коэффициент корреляции= .986

с параметром-(Фракция 0,01-0,005 мм в %)

коэффициент корреляции= .970

с параметром-(Фракция 0,005-0,001 мм в %)

коэффициент корреляции= .980

с параметром-(Фракция менее 0,001 мм в %)

коэффициент корреляции= .984

Результаты факторного анализа на ранжированной матрице исследования по 1963 году:

Таблица-17-Объединение по фактору 1.

в факторе- 1 базовый параметр- 6 (Фракция 0,05-0,01 мм в %)

Результаты факторного анализа на ранжированной матрице исследования по 1983 году:

Таблица-18-Объединение по фактору 1.

в факторе- 1 базовый параметр- 8

(Фракция 0,005-0,001 мм в %)

Согласно результатам корреляционного анализа для ранжированных матриц исследования имеем наличие сильных парных обусловленностей, потому что коэффициенты парной корреляции имеют значения по модулю больше 0,8.

Факторный анализ дал объединение всех параметров ранжированных нормализованных матриц исследования в одном факторе.

Это означает наличие одной групповой обусловленности для всех параметров исследования, одной природы их поведения.

Базовые параметры для ранжированных нормализованных матриц исследования:

-для матрицы исследования 1963 года базовый параметр-фракция 0,05-0,01 мм в %.

-для матрицы исследования 1983 года базовый параметр-фракция 0,005-0,001 мм в %.

Список литературы

1. Бендат Д. Ж., Пирсол А. Измерение и анализ случайных процессов. - М.: Мир, 2014.

2. Драйпер Н., Смит Г. Прикладной регрессионный анализ. - М.: Статистика, 1973.

3. Brandon D. B. Developing Mathematical Models for Computer Control, USA Journal, 2011, V.S,N7.

4. Харман Г. Современный факторный анализ.-М.:Сатистика, 1972.

5. Иберла К. Факторный анализ.-М.:Статистика, 1980.

6. Lawley D.M. The estimation of factor loadings by the method of maximum likelihood. Proc. roy. Soc. Edinb. Abo. 64-82(2010).

7. Kaiser H. F. [1]. The varimax criterio for analytic rotation in factor analysis. Psychometrica, 23, 187-200(1958).

Приложение

тексты программ

Набор FACT.CPP

//c ТЕКСТ ПРОГРАММЫ

//c ПО ФАКТОРНОМУ АНАЛИЗУ

//c (FACT.cpp)

#include<windows.h>

#include<string.h>

#include<stdio.h>

#include"koregmnk.h"

#include<math.h>

#include<alloc.h>

void factory()

{

extern void faca(long int n,long int m,float *x,long int *mnpi,long int ki);

extern char **anaz;

extern long int knaz;

extern FILE *fp6,*fp4,*fp5;

char prob=' ',pea;

long int n,m,ki,*mnpi,nm1,nm,i,j;

float *x;

fp6=fopen("factwiwc","w");

fp4=fopen("sled","r");

fp5=fopen("name","r");

anaz=(char **)farmalloc(501*4);

if(anaz==NULL)

{

fprintf(fp6,"Not memory for anaz i-vector\n");

return;

}

for(j=0;j<=500;j++)

{

*(anaz+j)=(char *)farmalloc(61*1);

if(*(anaz+j)==NULL)

{

fprintf(fp6,"Not memory for anaz-vector\n");

return;

}

}

x=(float *)farmalloc(250001*4);

if(x==NULL)

{

fprintf(fp6,"Not memory for x vector\n");

return;

}

mnpi=(long int *)farmalloc(501*4);

if(mnpi==NULL)

{

fprintf(fp6,"Not memory for mnpi vector\n");

return;

}

fscanf(fp5,"%ld",&knaz);

//fprintf(fp6,"%ld\n",knaz);

fgetc(fp5);/*чтение \n*/

//fprintf(fp6,"%ld\n",knaz);

if(knaz<=500) goto m6721;

fprintf(fp6,"количество параметров-столбиков >500\n");

return;

m6721:;

for(i=1;i<=knaz;i++)

for(j=1;j<=60;j++)

anaz[i][j]=' ';

for(i=1;i<=knaz;i++)

{

for(j=1;j<=61;j++)

{

pea=fgetc(fp5);

//fprintf(fp6,"i=%ld j=%ld pea=%c\n",i,j,pea);

if(pea=='\n') goto man;

anaz[i][j]=pea;

// fprintf(fp6,"%c",anaz[i][j]);

}

man:;

// fprintf(fp6,"\n");

}

for(i=1;i<=81;i++)

{

pea=fgetc(fp4);

if(pea=='\n') goto conti;

fprintf(fp6,"%c",pea);

}

conti:;

fprintf(fp6,"\n");

fscanf(fp4,"%ld%ld",&n,&m);

//fprintf(fp6,"n=%ld m=%ld\n",n,m);

nm=n*m;

if(nm<=250000) goto m6781;

fprintf(fp6,"размер матрицы больше 250000\n");

return;

m6781:;

for(i=1;i<=n;i++)

{

for(j=i;j<=nm;j=j+n)

{

fscanf(fp4,"%f",&x[j]);

//fprintf(fp6,"%f ",x[j]);

}

// fprintf(fp6,"\n");

}

fprintf(fp6,"\n");

for(j=1;j<=m;j++)

mnpi[j]=j;

ki=m;

faca(n,m,x,mnpi,ki);

fclose(fp5);

fclose(fp4);

fclose(fp6);

/*farfree(anaz);

farfree(mnpi);

farfree(x); */

return;

}

void faca(long int n,long int m,float *x,long int *mnpi,long int kis)

{

extern void corre(long int n,long int m,long int id,

float *x,float *xbar,float *std,

float *v,float *r,float *b,float *d,float *t);

extern void eigen(float *r,float *rd,long int m,long int mv);

extern void trace(long int m,float *r,float con,long int *pk,float *d);

extern void load(long int m,long int k1,float *r,float *v);

extern void gf(long int m,long int k1,float *v);

extern void grpk(long int m,long int k1,float *v,long int *mnpy,

long int ip1,long int *mnpi,long int kis);

extern void varmx(long int m,long int k1,float *v,long int *pnc,float *tv,float *h,

float *f,float *d);

extern char **anaz;

extern long int knaz;

extern FILE *fp6;

long int *mnpy,*mnp,*isumf,*inomf,iwkf=1,kj3,nz,iz,il,k,ipe,nc,

k1,id,i,mv,km,mm,ird,i1,n1,n2,mtot,ipo,ip1,inab,jpo,

mmk,ifac,in,ik,j,mst,kk1,is,nm,mn1,nm1,ij,ijk,mne,mpe,mke;

float *zmaxf,*std,*b,*xbar,*t,*d,*tv,*h,*f,*bsr,*v,

*r,*rd,con,z,xmax;

r=(float *)farmalloc(250001*4);

if(r==NULL)

{

fprintf(fp6,"Not memory for r vector\n");

return;

}

v=(float *)farmalloc(250001*4);

if(v==NULL)

{

fprintf(fp6,"Not memory for v vector\n");

return;

}

rd=(float *)farmalloc(250001*4);

if(rd==NULL)

{

fprintf(fp6,"Not memory for rd vector\n");

return;

}

mnpy=(long int *)farmalloc(501*4);

if(mnpy==NULL)

{

fprintf(fp6,"Not memory for mnpy vector\n");

return;

}

mnp=(long int *)farmalloc(501*4);

if(mnp==NULL)

{

fprintf(fp6,"Not memory for mnp vector\n");

return;

}

isumf=(long int *)farmalloc(501*4);

if(isumf==NULL)

{

fprintf(fp6,"Not memory for isumf vector\n");

return;

}

inomf=(long int *)farmalloc(501*4);

if(inomf==NULL)

{

fprintf(fp6,"Not memory for inomf vector\n");

return;

}

zmaxf=(float *)farmalloc(501*4);

if(zmaxf==NULL)

{

fprintf(fp6,"Not memory for zmaxf vector\n");

return;

}

std=(float *)farmalloc(501*4);

if(std==NULL)

{

fprintf(fp6,"Not memory for std vector\n");

return;

}

b=(float *)farmalloc(501*4);

if(b==NULL)

{

fprintf(fp6,"Not memory for b vector\n");

return;

}

xbar=(float *)farmalloc(501*4);

if(xbar==NULL)

{

fprintf(fp6,"Not memory for xbar vector\n");

return;

}

t=(float *)farmalloc(501*4);

if(t==NULL)

{

fprintf(fp6,"Not memory for t vector\n");

return;

}

d=(float *)farmalloc(501*4);

if(d==NULL)

{

fprintf(fp6,"Not memory for d vector\n");

return;

}

tv=(float *)farmalloc(501*4);

if(tv==NULL)

{

fprintf(fp6,"Not memory for tv vector\n");

return;

}

h=(float *)farmalloc(501*4);

if(h==NULL)

{

fprintf(fp6,"Not memory for h vector\n");

return;

}

f=(float *)farmalloc(501*4);

if(f==NULL)

{

fprintf(fp6,"Not memory for f vector\n");

return;

}

bsr=(float *)farmalloc(501*4);

if(bsr==NULL)

{

fprintf(fp6,"Not memory for bsr vector\n");

return;

}

con=0.6;

kj3=m;

nz=1;

iz=m*(m+1)/2;

//fprintf(fp6,"\ndo m=%ld iz=%ld\n",m,iz);

il=0;

k1=3;

id=1;

corre(n,m,id,x,xbar,std,v,r,b,d,t);

//fprintf(fp6,"\nr m=%ld iz=%ld\n",m,iz);

//for(i=1;i<=iz;i++)

//fprintf(fp6,"%7.3f ",r[i]);

//fprintf(fp6,"\n");

for(i=1;i<=m;i++)

{

bsr[i]=xbar[i];

b[i]=xbar[i];

}

mv=0;

eigen(r,rd,m,mv);

//fprintf(fp6,"\np ei r\n");

//for(i=1;i<=iz;i++)

//fprintf(fp6,"%7.3f ",r[i]);

//fprintf(fp6,"\n rd\n");

mm=m*m;

//for(i=1;i<=mm;i++)

//fprintf(fp6,"%f\n",rd[i]);

il=0;

id=1;

km=1;

m7007: il=il+1;

xbar[il]=r[km];

id=id+1;

km=km+id;

if(km<=iz) goto m7007;

il=0;

mm=m*m;

ird=0;

for(jpo=1;jpo<=m;jpo++)

{

//fprintf(fp6,"\n собственное значение=%12.4f\n",xbar[jpo]);

il=0;

for(ipo=jpo;ipo<=mm;ipo=ipo+m)

{

il=il+1;

h[il]=rd[ipo];

//fprintf(fp6,"\n собственный вектор(столбик)\n");

//for(ipo=1;ipo<=m;ipo++)

//fprintf(fp6,"%f\n",h[ipo]);

}}

n1=1;

n2=m;

il=0;

trace(m,r,con,&k,d);

k1=k;

iz=m*m;

mtot=m;

for(il=1;il<=iz;il++)

v[il]=rd[il];

load(m,k1,r,v);

m=mtot;

for(il=1;il<=iz;il++)

rd[il]=v[il];

if(iwkf==0) goto m70280;

mm=m*m;

for(ipo=1;ipo<=m;ipo++)

{

il=0;

for(jpo=ipo;jpo<=mm;jpo=jpo+m)

{

il=il+1;

h[il]=v[jpo];

}

}

fprintf(fp6,"\nдо вращения число факторов-%ld\n",k1);

gf(m,k1,v);

for(i=1;i<=m;i++)

mnpy[i]=0;

ip1=0;

grpk(m,k1,v,mnpy,ip1,mnpi,kis);

fprintf(fp6,"\n Попадание наблюдений в факторы по максимальным\n");

fprintf(fp6,"значениям произведений факторных нагрузок на значения\n");

fprintf(fp6,"факторов(до вращения) из всех параметров в наблюдении\n");

inab=0;

m=kj3;

mm=m*m;

for(j=1;j<=m;j++)

isumf[j]=0;

ik=1;

m7027: iz=ik;

il=0;

m7024: il=il+1;

if(il>kj3) goto m7023;

m7022:;

t[il]=x[iz];

iz=iz+n;

goto m7024;

m7023:;

//c-вычисление значений факторов

il=0;

for(jpo=1;jpo<=k1;jpo++)

{

d[jpo]=0;

for(ipo=1;ipo<=m;ipo++)

{

il=il+1;

d[jpo]=d[jpo]+v[il]*t[ipo];

}

d[jpo]=d[jpo]/xbar[jpo];

}

mtot=m;

inab=inab+1;

//c-проверка

mmk=m*k1;

for(jpo=1;jpo<=m;jpo++)

{

il=0;

f[jpo]=0;

for(ipo=jpo;ipo<=mmk;ipo=ipo+m)

{

il=il+1;

z=d[il];

h[il]=v[ipo]*z;

rd[ipo]=h[il];

f[jpo]=f[jpo]+ v[ipo]*z;

}

xmax=fabs(h[1]);

ifac=1;

for(ipo=1;ipo<=k1;ipo++)

{

if(xmax>=fabs(h[ipo])) goto m8621;

xmax=fabs(h[ipo]);

ifac=ipo;

m8621:;

}

zmaxf[jpo]=xmax;

inomf[jpo]=ifac;

}

xmax=fabs(zmaxf[1]);

in=inomf[1];

for(jpo=1;jpo<=m;jpo++)

{

if(xmax>=fabs(zmaxf[jpo])) goto m33336;

xmax=fabs(zmaxf[jpo]);

in=inomf[jpo];

m33336:;

}

isumf[in]=isumf[in]+1;

fprintf(fp6," в факторе-%ld наблюдение-%ld\n",in,inab);

for(i=1;i<=m;i++)

mnpy[i]=0;

ip1=0;

ik=ik+1;

if(ik<=n) goto m7027;

fprintf(fp6,"\n Распределение наблюдений по факторам\n");

for(j=1;j<=k1;j++)

{

if(isumf[j]==0) goto m33333;

fprintf(fp6," фактор-%ld ,число наблюдений в этом факторе=%ld\n",

j,isumf[j]);

m33333:;

}

m70280:;

fprintf(fp6,"\nпосле вращения\n");

//c-Нахождение номеров нулевых строчек в матрице факторных нагрузок v()

//c-массивы MNP(),MNPY()-номера нулевых строчек в матрице v(),верно!

mst=m;

kk1=k1*m;

ip1=0;

for(i=1;i<=m;i++)

{

is=0;

for(j=i;j<=kk1;j=j+m)

{

is=is+1;

b[is]=v[j];

}

for(j=1;j<=is;j++)

if(b[j]!=0) goto m101;

ip1=ip1+1;

mnp[ip1]=i;

mnpy[ip1]=i;

m101:;

}

if(ip1==0) goto m104;

fprintf(fp6,"\n otc\n");

for(i=1;i<=ip1;i++)

fprintf(fp6,"%ld\n",mnp[i]);

//c-Уплотнение по нулевым строчкам в матрице факторных нагрузок v()

//c-элементы в матрице идут в массиве v() по столбцам,верно!

for(i=1;i<=ip1;i++)

{

nm=mnp[i];

if(nm==m) goto m106;

mn1=nm+1;

nm1=mn1;

for(j=nm1;j<=m;j=j+1)

{

is=0;

for(ij=j;ij<=kk1;ij=ij+m)

{

is=is+1;

b[is]=v[ij];

}

is=0;

ij=j-1;

for(ijk=ij;ijk<=kk1;ijk=ijk+m)

{

is=is+1;

v[ijk]=b[is];

}

}

for(ij=i;ij<=ip1;ij++)

mnp[ij]=mnp[ij]-1;

}

m=m-ip1;

m104:;

goto m118;

m106:;

m=m-1;

m118:;

//c-формирование матрицы факторных нагрузок после уплотнения строк

//c-элементы матрицы расположены по столбцам,верно!

for(i=2;i<=k1;i++)

{

mne=(i-1)*m;

mpe=(i-1)*mst+1;

mke=mpe+m-1;

for(ipe=mpe;ipe<=mke;ipe++)

{

mne=mne+1;

v[mne]=v[ipe];

}}

varmx(m,k1,v,&nc,tv,h,f,d);

fprintf(fp6,"число факторов=%ld\n",k1);

gf(m,k1,v);

grpk(m,k1,v,mnpy,ip1,mnpi,kis);

/*farfree(r);

farfree(v);

farfree(rd);

farfree(mnpy);

farfree(mnp);

farfree(isumf);

farfree(inomf);

farfree(zmaxf);

farfree(std);

farfree(b);

farfree(xbar);

farfree(t);

farfree(d);

farfree(bsr);

farfree(tv);

farfree(h);

farfree(f); */

return;

}

void gf(long int m,long int k1,float *v)

{

extern FILE *fp6;

longint md[6],jk,in,ir,i,j,l;

float s1[501],s2[501];

jk=k1*m;

in=1;

ir=m;

for(i=1;i<=5;i++)

md[i]=i;

for(j=1;j<=k1;j++)

s2[j]=0.;

for(i=1;i<=m;i++)

{

l=0;

for(j=i;j<=jk;j=j+m)

{

l=l+1;

s1[l]=v[j]*v[j];

s2[l]=s2[l]+v[j]*v[j];

}

}

fprintf(fp6,"\n Таблица-Сумма квадратов нагрузок по факторам\n");

for(i=1;i<=49;i++)

{

if(i<=9) fprintf(fp6," ");

else

fprintf(fp6,"-");

}

fprintf(fp6,"\n |НОМЕР ФАКТОРА|СУММА КВАДРАТОВ НАГРУЗОК|\n");

for(i=1;i<=49;i++)

{

if(i<=9) fprintf(fp6," ");

else

fprintf(fp6,"-");

}

fprintf(fp6,"\n");

for(j=1;j<=k1;j++)

{

fprintf(fp6," |%13d|%24.3f|\n",j,s2[j]);

for(i=1;i<=49;i++)

{

if(i<=9) fprintf(fp6," ");

else

fprintf(fp6,"-");

}

fprintf(fp6,"\n");

}

return;

}

void corre(long int pn,long int pm,long int pid,float *x,float *xbar,float *std,

float *rx,float *r,float *b,float *d,float *t)

{

long int i,j,kk,k,l,jk,io,n,m;

float fn,fkk,mmm;

n=pn;

m=pm;

io=pid;

for(j=1;j<=m;j++)

{

b[j]=0.0;

t[j]=0.0;

}

k=(m*m+m)/2;

for(i=1;i<=k;i++)

r[i]=0.0;

fn=n;

l=0;

//if(io!=0) goto m105; else goto m127;

m105:for(j=1;j<=m;j++)

{

for(i=1;i<=n;i++)

{

l++;

t[j]=t[j]+x[l];

}

xbar[j]=t[j];

t[j]=t[j]/fn;

}

for(i=1;i<=n;i++)

{

jk=0;

l=i-n;

for(j=1;j<=m;j++)

{

l=l+n;

d[j]=x[l]-t[j];

b[j]=b[j]+d[j];

}

for(j=1;j<=m;j++)

for(k=1;k<=j;k++)

{

jk++;

r[jk]=r[jk]+d[j]*d[k];

}

}

m205:jk=0;

for(j=1;j<=m;j++)

{

xbar[j]=xbar[j]/fn;

for(k=1;k<=j;k++)

{

jk++;

r[jk]=r[jk]-b[j]*b[k]/fn;

}

}

jk=0;

for(j=1;j<=m;j++)

{

jk=jk+j;

mmm=fabs(r[jk]);

std[j]=sqrt(mmm);

}

for(j=1;j<=m;j++)

for(k=j;k<=m;k++)

{

jk=j+(k*k-k)/2;

l=m*(j-1)+k;

rx[l]=r[jk];

l=m*(k-1)+j;

rx[l]=r[jk];

if((std[j]*std[k])!=0) goto m225; else goto m222;

m222:r[jk]=0.0;

goto m230;

m225:r[jk]=r[jk]/(std[j]*std[k]);

m230:;

}

fn=sqrt(fn-1);

for(j=1;j<=m;j++)

std[j]=std[j]/fn;

l=-m;

for(i=1;i<=m;i++)

{

l=l+m+1;

b[i]=rx[l];

}

//pn=n;

//pm=m;

//pid=io;

return;

}

void eigen(float *a,float *r,long int n,long int mv)

{

longint iq,j,ij,i,ia,ind,l,m,mq,lq,lm,ll,mm,ilq,jq,k,

imq,im,il,ilr,imr,jk;

float range,anorm,anrmx,thr,x,y,sinx,sinx2,cosx,cosx2,sincs;

m5: range=1.0e-6;

if(mv!=1) goto m10;

else

goto m25;

m10: iq=-n;

for(j=1;j<=n;j++)

{

iq=iq+n;

for(i=1;i<=n;i++)

{

ij=iq+i;

r[ij]=0.0;

if(i!=j) goto m20;

else

goto m15;

m15: r[ij]=1.0;

m20:;

}

}

m25: anorm=0.0;

for(i=1;i<=n;i++)

for(j=i;j<=n;j++)

{

if(i!=j) goto m30;

else

goto m35;

m30: ia=i+(j*j-j)/2;

anorm=anorm+a[ia]*a[ia];

m35:;

}

if(anorm<=0) goto m165;

else

goto m40;

m40: anorm=1.414*pow(anorm,0.5);

anrmx=anorm*range/(1.*n);

ind=0;

thr=anorm;

m45: thr=thr/(1.*n);

m50: l=1;

m55: m=l+1;

m60: mq=(m*m-m)/2;

lq=(l*l-l)/2;

lm=l+mq;

m62: if(fabs(a[lm])>=thr) goto m65;

else

goto m130;

m65: ind=1;

ll=l+lq;

mm=m+mq;

x=0.5*(a[ll]-a[mm]);

m68: y=-a[lm]/pow(a[lm]*a[lm]+x*x,0.5);

if(x>=0) goto m75;

else

goto m70;

m70: y=-y;

m75: sinx=y/pow(2.*(1.+(pow(1.-y*y,0.5))),0.5);

sinx2=sinx*sinx;

m78: cosx=pow(1.0-sinx2,0.5);

cosx2=cosx*cosx;

sincs=sinx*cosx;

ilq=n*(l-1);

imq=n*(m-1);

for(i=1;i<=n;i++)

{

iq=(i*i-i)/2;

if(i!=l) goto m80;

else

goto m115;

m80: if(i<m) goto m85;

if(i==m) goto m115;

if(i>m) goto m90;

m85: im=i+mq;

goto m95;

m90: im=m+iq;

m95: if(i>=l) goto m105;

else

goto m100;

m100: il=i+lq;

goto m110;

m105: il=l+iq;

m110: x=a[il]*cosx-a[im]*sinx;

a[im]=a[il]*sinx+a[im]*cosx;

a[il]=x;

m115:if(mv!=1) goto m120;

else

goto m125;

m120: ilr=ilq+i;

imr=imq+i;

x=r[ilr]*cosx-r[imr]*sinx;

r[imr]=r[ilr]*sinx+r[imr]*cosx;

r[ilr]=x;

m125:;

}

x=2.0*a[lm]*sincs;

y=a[ll]*cosx2+a[mm]*sinx2-x;

x=a[ll]*sinx2+a[mm]*cosx2+x;

a[lm]=(a[ll]-a[mm])*sincs+a[lm]*(cosx2-sinx2);

a[ll]=y;

a[mm]=x;

m130: if(m!=n) goto m135;

else

goto m140;

m135: m=m+1;

goto m60;

m140:if(l!=(n-1)) goto m145;

else

goto m150;

m145: l=l+1;

goto m55;

m150: if(ind!=1) goto m160;

else

goto m155;

m155: ind=0;

goto m50;

m160: if(thr<=anrmx) goto m165;

else

goto m45;

m165: iq=-n;

for(i=1;i<=n;i++)

{

iq=iq+n;

ll=i+(i*i-i)/2;

jq=n*(i-2);

for(j=i;j<=n;j++)

{

jq=jq+n;

mm=j+(j*j-j)/2;

if(a[ll]>=a[mm]) goto m185;

else

goto m170;

m170: x=a[ll];

a[ll]=a[mm];

a[mm]=x;

if(mv!=1) goto m175;

else

goto m185;

m175: for(k=1;k<=n;k++)

{

ilr=iq+k;

imr=jq+k;

x=r[ilr];

r[ilr]=r[imr];

m180: r[imr]=x;

}

m185:;

}}

return;

}

void load(long int m,long int k,float *r,float *v)

{

longint l,jj,j,i;

float sq;

l=0;

jj=0;

for(j=1;j<=k;j++)

{

jj=jj+j;

m150: sq=pow(r[jj],0.5);

for(i=1;i<=m;i++)

{

l=l+1;

m160: v[l]=sq*v[l];

}

}

return;

}

void trace(long int m,float *r,float con,long int *pk,

float *d)

{

longint i,l,k;

float fm;

fm=m;

l=0;

for(i=1;i<=m;i++)

{

l=l+i;

m100: d[i]=r[l];

}

k=0;

for(i=1;i<=m;i++)

{

if(d[i]>=con) goto m105;

else

goto m120;

m105: k=k+1;

m110: d[i]=d[i]/fm;

}

m120: for(i=2;i<=k;i++)

m130: d[i]=d[i]+d[i-1];

*pk=k;

return;

}

void varmx(long int m,long int k,float *a,long int *pnc,

float *tv,float *h,float *f,

float *d)

{

longint ll,nv,nc,i,j,l,lb,l1,ii,l2,l3,l4,k1;

float eps,tvlt,fn,ffn,cons,aa,bb,cc,dd,u,t,b,cos4t,

sin4t,tan4t,sinp,cosp,ctn4t,cos2t,sin2t,cost,sint;

eps=0.00116;

tvlt=0.0;

ll=k-1;

nv=1;

nc=0;

fn=m;

ffn=fn*fn;

cons=0.7071066;

for(i=1;i<=m;i++)

{

h[i]=0.0;

for(j=1;j<=k;j++)

{

l=m*(j-1)+i;

m110: h[i]=h[i]+a[l]*a[l];

}}

for(i=1;i<=m;i++)

{

m115: h[i]=pow(h[i],0.5);

for(j=1;j<=k;j++)

{

l=m*(j-1)+i;

m120: a[l]=a[l]/h[i];

}

}

goto m132;

m130: nv=nv+1;

tvlt=tv[nv-1];

m132: tv[nv]=0.0;

for(j=1;j<=k;j++)

{

aa=0.0;

bb=0.0;

lb=m*(j-1);

for(i=1;i<=m;i++)

{

l=lb+i;

cc=a[l]*a[l];

aa=aa+cc;

m140: bb=bb+cc*cc;

}

m150: tv[nv]=tv[nv]+(fn*bb-aa*aa)/ffn;

}

if(nv>=51) goto m430;

m160: if((tv[nv]-tvlt)<=(1.e-7)) goto m170;

else

goto m190;

m170: nc=nc+1;

if(nc<=3) goto m190;

else goto m430;

m190: for(j=1;j<=ll;j++)

{

l1=m*(j-1);

ii=j+1;

for(k1=ii;k1<=k;k1++)

{

l2=m*(k1-1);

aa=0.0;

bb=0.0;

cc=0.0;

dd=0.0;

for(i=1;i<=m;i++)

{

l3=l1+i;

l4=l2+i;

u=(a[l3]+a[l4])*(a[l3]-a[l4]);

t=a[l3]*a[l4];

t=t+t;

cc=cc+(u+t)*(u-t);

dd=dd+2.0*u*t;

aa=aa+u;

m230: bb=bb+t;

}

t=dd-2.0*aa*bb/fn;

b=cc-(aa*aa-bb*bb)/fn;

if(t<b) goto m280;

if(t==b) goto m240;

if(t>b) goto m320;

m240: if((t+b)<eps) goto m420;

m250: cos4t=cons;

sin4t=cons;

goto m350;

m280: tan4t=fabs(t)/ fabs(b);

if(tan4t<eps) goto m300;

m290: cos4t=1.0/ pow(1.0+tan4t*tan4t,0.5);

sin4t=tan4t*cos4t;

goto m350;

m300:if(b>=0) goto m420;

m310: sinp=cons;

cosp=cons;

goto m400;

m320: ctn4t= fabs(t/b);

if(ctn4t<eps) goto m340;

m330: sin4t=1.0/ pow(1.0+ctn4t*ctn4t,0.5);

cos4t=ctn4t*sin4t;

goto m350;

m340: cos4t=0.0;

sin4t=1.0;

m350: cos2t=pow((1.0+cos4t)/2.0,0.5);

sin2t=sin4t/(2.0*cos2t);

m355: cost=pow((1.0+cos2t)/2.0,0.5);

sint=sin2t/(2.0*cost);

if(b<=0) goto m370;

m360: cosp=cost;

sinp=sint;

goto m380;

m370: cosp=cons*cost+cons*sint;

m375: sinp=fabs(cons*cost-cons*sint);

m380: if(t>0) goto m400;

m390: sinp=-sinp;

m400: for(i=1;i<=m;i++)

{

l3=l1+i;

l4=l2+i;

aa=a[l3]*cosp+a[l4]*sinp;

a[l4]=-a[l3]*sinp+a[l4]*cosp;

m410: a[l3]=aa;

}

m420:;

}}

goto m130;

m430:for(i=1;i<=m;i++)

{

for(j=1;j<=k;j++)

{

l=m*(j-1)+i;

m440: a[l]=a[l]*h[i];}}

nc=nv-1;

for(i=1;i<=m;i++)

m450: h[i]=h[i]*h[i];

for(i=1;i<=m;i++)

{

f[i]=0.0;

for(j=1;j<=k;j++)

{

l=m*(j-1)+i;

m460: f[i]=f[i]+a[l]*a[l];

}

m470: d[i]=h[i]-f[i];

}

*pnc=nc;

return;

}

void grpk(long int m,long int k1,float *v,long int *mnp,long int ip1,

long int *mnpi,long int kis)

{

extern FILE *fp6;

extern char **anaz;

extern long int knaz;

longint mnf[501],mrez[501],nf[501],mgp[501],iti,

ikr,i,is,ip,kt,ins,jst,nft,i1,ib,j,ir,ismax,it,ii;

float zna[501],xmax,s1;

char prob=' ';

ikr=0;

if(ip1!=0) goto m91;

for(i=1;i<=kis;i++)

m92: mrez[i]=mnpi[i];

goto m93;

m91:;

for(is=1;is<=kis;is++)

{

for(ip=1;ip<=ip1;ip++)

if(is==mnp[ip]) goto m1;

ikr=ikr+1;

mrez[ikr]=mnpi[is];

m1:;

}

m93:;

kt=m*k1;

for(i=1;i<=m;i++)

{

xmax=fabs(v[i]);

ins=i;

for(j=i;j<=kt;j=j+m)

{

if(xmax>=fabs(v[j])) goto m6;

xmax=fabs(v[j]);

ins=j;

m6:;

}

jst=ins/m+1;

s1=ins*1./m;

i1=s1;

if(s1==i1) jst=i1;

zna[i]=v[ins];

nf[i]=jst;

m5:;

}

for(i=1;i<=m;i++)

{

nft=nf[i];

mnf[i]=nft;

if(i==1) goto m201;

ii=i-1;

for(ir=1;ir<=ii;ir++)

if(nft==mnf[ir]) goto m10;

m201:;

ib=0;

for(j=i;j<=m;j++)

{

if(nft!=nf[j]) goto m11;

ib=ib+1;

mgp[ib]=j;

m11:;}

fprintf(fp6,"\nТаблица-Объединение по фактору-%ld\n",nft);

for(iti=1;iti<=58;iti++)

fprintf(fp6,"-");

fprintf(fp6,"\n");

fprintf(fp6,"|НОМЕР| НАЗВАНИЕ ПАРАМЕТРА |");

fprintf(fp6," НАГРУЗКА |\n");

for(iti=1;iti<=58;iti++)

fprintf(fp6,"-");

fprintf(fp6,"\n");

ir=mgp[1];

ismax=mrez[ir];

xmax=fabs(zna[ir]);

for(it=1;it<=ib;it++)

{

ir=mgp[it];

is=mrez[ir];

if(xmax>=fabs(zna[ir])) goto m12102;

ismax=is;

xmax=fabs(zna[ir]);

m12102:;

for(j=35;j<=60;j++)

{

if(anaz[is][j]==prob) goto m214;

goto m215;

m214:;}

fprintf(fp6,"|%5d|",is);

for(j=1;j<=34;j++)

fprintf(fp6,"%c",anaz[is][j]);

fprintf(fp6,"|%11.4f |\n",zna[ir]);

for(iti=1;iti<=58;iti++)

fprintf(fp6,"-");

fprintf(fp6,"\n");

goto m13;

m215:;

fprintf(fp6,"|%5d|",is);

for(j=1;j<=34;j++)

fprintf(fp6,"%c",anaz[is][j]);

fprintf(fp6,"| |\n");

fprintf(fp6,"| |");

for(j=35;j<=60;j++)

fprintf(fp6,"%c",anaz[is][j]);

fprintf(fp6," |%11.4f |\n",zna[ir]);

for(iti=1;iti<=58;iti++)

fprintf(fp6,"-");

fprintf(fp6,"\n");

m13:;}

fprintf(fp6,"в факторе-%5d базовый параметр-%5d\n",nft,ismax);

for(j=1;j<=60;j++)

fprintf(fp6,"%c",anaz[ismax][j]);

fprintf(fp6,"\n");

m10:;}

return;}

Набор menu.cpp(главная и оконная функции при нормализации)

/*WORK WITH MENU NO ACCELERATOR ¦L++TL T ¦+=¦ ++¦ L¦T+T+¦LT+¦+T*/

#include<windows.h>

#include<string.h>

#include<stdio.h>

#include<alloc.h>

#include<dos.h>

#include"menu.h"

FILE *fp5,*fp6,*fp7,*fp8,*fp9,*fpo,*fp4,*fp3;

LRESULT CALLBACK WindowFunc(HWND,UINT,WPARAM,LPARAM);

char szWinName[]="MY WINDOW ¦++ юъэю";/*NAME OF CLASS OF WINDOWS шь ъырёёр юъэр*/

int WINAPI WinMain(HINSTANCE hThisInst,

HINSTANCE hPrevInst,

LPSTR lpszArgs,

int nWinMode)

{

HWND hwnd;

MSG msg;

WNDCLASS wcl;

/*DEFINE CLASS OF WINDOWSа*/

wcl.hInstance=hThisInst;/*DESCRIPTOR OF APPLICATIONS фхёъЁшяЄюЁ яЁшыюцхэш */

wcl.lpszClassName=szWinName;/* NAME OF CLASSES OF WINDOWS шь ъырёёр юъэр*/

wcl.lpfnWndProc=WindowFunc;/*FUNCTION OF WINDOWS ЇєэъЎш юъэр*/

wcl.style=0; /*STYLE ёЄшы№ яю єьюыўрэш¦*/

wcl.hIcon=LoadIcon(NULL,IDI_APPLICATION);/*STANDART ICON ёЄрэфрЁЄэр шъюэър*/

wcl.hCursor=LoadCursor(NULL,IDC_ARROW);/*STANDART CURSOR ёЄрэфрЁЄэvщ ъєЁёюЁ*/

wcl.lpszMenuName="MYMENU";/*NAME OF MENUS шь ьхэ¦*/

wcl.cbClsExtra=0;/*NOT INFORMATION эхЄ фюяюыэшЄхы№эющ */

wcl.cbWndExtra=0;/*NOT INFORMATION шэЇюЁьрЎшш*/

wcl.hbrBackground=(HBRUSH)GetStockObject(WHITE_BRUSH);

/* REGISTRATIVE CLASS OF WINDOWS ЁхушёЄЁрЎш ъырёёр юъэр*/

if(!RegisterClass(&wcl))

return 0;

/* CREATE WINDOW */

hwnd=CreateWindow(szWinName,/*NAME OF CLASS OF WINDOWS шь ъырёёр юъэр*/

"postroenie norm matric build6",/* CHAPTER чруюыютюъ*/

WS_OVERLAPPEDWINDOW,/* STYLE OF WINDOWS TTLT- +¦=L*/

CW_USEDEFAULT,/* X-COORDINATE DEFINE Windows ¦++¦-L=LTL +¦¦+-+T-+TT- WINDOWS*/

CW_USEDEFAULT,/* Y-COORDINATE DEFINE Windows*/

CW_USEDEFAULT,/* WIGTH DEFINE Windows +L¦L=L +¦¦+-+T-+TT- WINDOWS*/

CW_USEDEFAULT,/* HEGHT DEFINE WINDOWS тvёюЄр юяЁхфхы хЄё Windows*/

HWND_DESKTOP,/*NO PAORENT эхЄ ЁюфшЄхы№ёъюую юъэр */

NULL,/* NOT MENU эхЄ ьхэ¦*/

hThisInst,/*DESCRIPTOR OF APPLICATION фхёъЁшяЄюЁ яЁшыюцхэш */

NULL/* NOT ARGYMENTS схч фюяюыэшЄхы№эvї рЁуєьхэЄют*/

);

/*SHOW WINDOW яюърчрЄ№ юъэю ш яхЁхЁшёютрЄ№ ёюфхЁцшьюх */

ShowWindow(hwnd,nWinMode);

UpdateWindow(hwnd);

/* START CYCLE чряєёЄшЄ№ Ўшъы юсЁрсюЄъш ёююс•хэшщ*/

while(GetMessage(&msg,NULL,0,0))

{

TranslateMessage(&msg);/*KEYBORD ЁрчЁх°шЄ№ шёяюы№чютрэшх ъыртшрЄєЁv*/

DispatchMessage(&msg);/*TO WINDOWS тхЁэєЄ№ єяЁртыхэшх Windows*/

}

return msg.wParam;

}

/* FUNCTION OF WINDOWS ЇєэъЎш юъэр */

LRESULT CALLBACK WindowFunc(HWND hwnd,UINT message,

WPARAM wParam,

LPARAM lParam)

{

extern FILE *fpo;

extern long int soch();

int response;

long int ikone;

//fpo=fopen("okno","w");

switch(message)

{

case WM_COMMAND:

switch(LOWORD(wParam))

{

case IDM_RUN:

//fpo=fopen("okno","w");

//fprintf(fpo,"\n idrun pered soch \n");

MessageBox(hwnd," START","RUN",MB_OK);

//fpo=fopen("okno","w");

//fprintf(fpo,"\npered soch\n");

ikone=soch();

//fprintf(fpo,"\n posle soch ikone=%ld\n",ikone);

//fclose(fpo);

if(ikone==1)

MessageBox(hwnd," END RUN ","FINISH ",MB_OK);

else

MessageBox(hwnd," CONTINUE RUN ","NOT END ",MB_OK);

break;

}

break;

case WM_DESTROY:/*END OF PROGRAMMS-чртхЁ°хэшх яЁюуЁрььv*/

PostQuitMessage(0);

break;

default:

/*ALL MESSAGES TO WINDOWS тёх ёююс•хэш ,эх юсЁрсрЄvтрхьvх т фрээющ ЇєэъЎшш,

эряЁрты ¦Єё эр юсЁрсюЄъє яю єьюыўрэш¦*/

return DefWindowProc(hwnd,message,wParam,lParam);

}

return 0;

}

c ТЕКСТ ПРОГРАММЫ

c ПО ФАКТОРНОМУ АНАЛИЗУ

c (FACTMAIN.FOR)

$large

dimension x(22801),mnpi(151),zgl(80)

common anaz(151,60),knaz

open(unit=6,file='factwiw',status='Unknown')

open(unit=4,file='sled',status='old')

open(unit=5,file='name',status='old')

data prob/' '/

read(5,*) knaz

cif(knaz.gt.70) goto 781

do 782 i=1,knaz

do 456 j=1,60

anaz(i,j)=prob

456continue

read(5,785)(anaz(i,j),j=1,60)

785format(60a1)

782continue

781continue

read(4,1700) zgl

1700format(80a1)

write(6,1701) zgl

1701format(1x,80a1)

read(4,*) n,m

nm=n*m

nm1=n*4

do 1 i=1,n

cwrite(6,21) i

21format(1x,'i=',i3)

read(4,*) (x(j),j=i,nm,n)

cwrite(*,22) (x(j),j=i,nm1,n)

c22format(1x,

c write(6,*)(x(j),j=i,nm,n)

1 continue

ki=m

do 3 j=1,m

3 mnpi(j)=j

cpause 'faca'

call faca(n,m,x,mnpi,ki)

close(unit=5,status='keep')

close(unit=4,status='keep')

close(unit=6,status='keep')

end

subroutine data

return

end

SUBROUTINE FACA(N,M,X,MNPI,KIS)

dimension x(1),isumf(200),zmaxf(151),inomf(151)

dimension std(151),mnp(151)

dimension b(151),xbar(151)

dimension t(151),d(151),tv(151)

dimension h(151),f(151)

dimension mnpy(151)

DIMENSION R(11476),V(22801),rd(22801)

DIMENSION MNPI(1),bsr(151)

common anaz(151,60),knaz

con=0.6

write(*,8761)

8761 format(1x,'если выводим вклады факторов,то ввести 1,нет-0')

read(*,*) iwkf

ciwkf=1

cwrite(6,66011) kis

66011format(1x,'faca wxod kis=',i3)

9292 FORMAT(F6.2)

21 FORMAT(2I4)

KJ3=M

NZ=1

22 FORMAT(F6.2,E8.2,4F6.2)

IZ=M*(M+1)/2IL=0

K1=3

ID=1

CALL CORRE(N,M,ID,X,XBAR,STD,V,R,B,D,T)

cwrite(6,7621)

7621format(1x,'r'/)

cwrite(6,7622) (r(i),i=1,iz)

7622format(1x,10f7.3/)

306 FORMAT(' ',6f10.5)

DO 7021 I=1,M

bsr(i)=xbar(i)

B(I)=XBAR(I)

7021 CONTINUE

7005 FORMAT(' ','средние')

7006 FORMAT(' ','sqrt(d)')

7003 FORMAT(' ','матрица r')

MV=0

CALL EIGEN(R,rd,M,MV)

cwrite(6,*) 'p eg r'

cwrite(6,*) (r(i),i=1,iz)

mm=m*m

cwrite(6,*) 'rd'

cwrite(6,*) (rd(i),i=1,mm)

IL=0

ID=1

KM=1

7007 IL=IL+1

XBAR(IL)=R(KM)

ID=ID+1

KM=KM+ID

IF(KM-IZ) 7007,7007,7008

7008 CONTINUE

7009 FORMAT(' ','собственние значения')

c o

c WRITE(6,7009)

c WRITE(6,306)(XBAR(I),I=1,IL)

il=0

mm=m*m

cwrite(6,4012)

ird=0

do 4013 jpo=1,m

cwrite(6,4014) xbar(jpo)

4014 format(1x,'собственное значение=',f12.4)

il=0

do 4015 ipo=jpo,mm,m

il=il+1

h(il)=rd(ipo)

4015 continue

cwrite(6,4012)

4012 format(1x,'собственные вектора(столбики)')

cwrite(6,306) (h(ipo),ipo=1,m)

4013continue

N1=1

N2=M

IL=0

CALL TRACE(M,R,CON,K,D)

7012 FORMAT(' ','накопленные отнош собств')

c WRITE(6,63871)K,m

63871 FORMAT(1X,'tr число собств>0.8',I5,' m=',i5)

c WRITE(6,7012)

c write(6,63872)(d(i),i=1,m)

63872 FORMAT(1X,5F10.3)

K1=k

IZ=M*M

mtot=m

do 203 il=1,iz

203 v(il)=rd(il)

CALL LOAD(M,K1,R,V)

cwrite(6,7821) m,k1

7821format(1x,'load m=',i5,' k1=',i5)

m=mtot

do 6071 il=1,iz

rd(il)=v(il)

cr(il)=v(il)

6071 continue

if(iwkf.eq.0) goto 70280

write(6,501)

cwrite(6,5012)

5012 format(/1x,'матрица факторных нагрузок(столбцы-факторы)'/)

mm=m*m

do 5013 ipo=1,m

il=0

do 5014 jpo=ipo,mm,m

il=il+1

h(il)=v(jpo)

5014 continue

cwrite(6,306) (h(il),il=1,m)

5013continue

9232 FORMAT(' число факторов',I5)

c WRITE(6,501)

write(6,9232) k1

CALL GF(M,K1,V)

do 6037 i=1,m

mnpy(i)=0

6037 continue

ip1=0

CALL GRPK(M,K1,V,MNPY,IP1,MNPI,KIS)

6028 CONTINUE

501 format(/' ','до вращения'/)

502 FORMAT(//' ','после вращения'//)

write(6,55551)

55551format(/1x,'Попадание наблюдений в факторы по максимальным'/1x,

*'значениям произведений факторных нагрузок на значения'/1x,

*'факторов(до вращения) из всех параметров в наблюдении'/)

inab=0

m=kj3

mm=m*m

do 33331 j=1,m

isumf(j)=0

33331continue

IK=1

7027 IZ=IK

IL=0

7024 IL=IL+1

IF(IL-KJ3) 7022,7022,7023

7022 CONTINUE

c if(std(il).eq.0) std(il)=0.0001

cT(IL)=(X(IZ)-Bsr(IL))/STD(IL)

T(IL)=(X(IZ))

c write(6,7201) t(il)

7201 format(1x,'t=',f15.4)

IZ=IZ+N

GOTO 7024

7023 CONTINUE

c-вычисление значений факторов

il=0

do 16601 jpo=1,k1

d(jpo)=0

do 6602 ipo=1,m

il=il+1

d(jpo)=d(jpo)+v(il)*t(ipo)

6602 continue

d(jpo)=d(jpo)/xbar(jpo)

16601continue

mtot=m

inab=inab+1

c-проверка

mmk=m*k1

do 8601 jpo=1,m

il=0

f(jpo)=0

do 8602 ipo=jpo,mmk,m

il=il+1

z=d(il)

h(il)=v(ipo)*z

rd(ipo)=h(il)

cwrite(6,9871) il,v(ipo),z,h(il)

9871format(1x,'фактор-',i3,'A=',f10.4,' F=',f10.4,' AF=',f10.4)

f(jpo)=f(jpo)+ v(ipo)*z

8602 continue

xmax=abs(h(1))

ifac=1

do 8621 ipo=1,k1

if(xmax.ge.abs(h(ipo))) goto 8621

xmax=abs(h(ipo))

ifac=ipo

8621 continue

cwrite(6,18603) inab

18603format(/1x,'в наблюдении-',i3)

cwrite(6,8604) (anaz(jpo,jpx),jpx=1,60)

8604 format(1x,'параметр-',60a1)

cwrite(6,8605) ifac

8605 format(1x,'имеет больший вклад от фактора-',i3/)

zmaxf(jpo)=xmax

inomf(jpo)=ifac

cstop

creturn

8603format(1x,'inab=',i3)

8601 continue

xmax=abs(zmaxf(1))

in=inomf(1)

do 33336 jpo=1,m

if(xmax.ge.abs(zmaxf(jpo))) goto 33336

xmax=abs(zmaxf(jpo))

in=inomf(jpo)

33336 continue

isumf(in)=isumf(in)+1

write(6,44441) in,inab

44441format(1x,'в факторе-',i5,' наблюдение-',i5)

cwrite(6,1341) inab

1341 format(/1x,'для наблюдения-',i3/1x,

*' объединение по max произведениям a(j,p)*f(p) по факторам'/)

do 16037 i=1,m

mnpy(i)=0

16037 continue

ip1=0

cCALL GRPK(M,K1,rd,MNPY,IP1,MNPI,KIS)

cwrite(6,8607) inab

do 8608 il=1,m

c write(6,8609) f(il),t(il)

8609 format(1x,'f=',f12.4,' t=',f12.4)

8608continue

8607format(1x,'prow inab=',i3/)

IK=IK+1

IF(IK-N) 7027,7027,7028

7028 CONTINUE

write(6,33332)

33332format(//1x,'Распределение наблюдений по факторам'/)

do 33333 j=1,k1

if(isumf(j).eq.0) goto 33333

write(6,33334) j,isumf(j)

33334format(1x,'фактор-',i3,' число наблюдений в этом факторе=',i5)

33333continue

70280continue

write(6,502)

c-Нахождение номеров нулевых строчек в матрице факторных нагрузок v()

c-массивы MNP(),MNPY()-номера нулевых строчек в матрице v(),верно!

MST=M

KK1=K1*M

IP1=0

DO 101 I=1,M

IS=0

DO 102 J=I,KK1,M

IS=IS+1

B(IS)=V(J)

102 CONTINUE

DO 103 J=1,IS

IF(B(J).NE.0) GOTO 101

103 CONTINUE

IP1=IP1+1

MNP(IP1)=I

MNPY(IP1)=I

101 CONTINUE

IF(IP1.EQ.0) GO TO 104

WRITE(6,115)(MNP(I),I=1,IP1)

115 FORMAT(1X,'OTC',10I5)

cwrite(6,66012) kis

66012format(1x,' faca prom kis=',i3)

c-Уплотнение по нулевым строчкам в матрице факторных нагрузок v()

c-элементы в матрице идут в массиве v() по столбцам,верно!

DO 105 I=1,IP1

NM=MNP(I)

IF(NM.EQ.M) GO TO 106

MN1=NM+1

NM1=MN1

DO 107 J=NM1,M

IS=0

DO 108 IJ=J,KK1,M

IS=IS+1

B(IS)=V(IJ)

108 CONTINUE

cwrite(6,66021) kis

66021format(1x,'faca posle 108 kis=',i3)

IS=0

IJ=J-1

DO 109 IJK=IJ,KK1,M

IS=IS+1

V(IJK)=B(IS)

109 CONTINUE

107 CONTINUE

cwrite(6,66022) kis

66022format(1x,'faca posle 107 kis=',i3)

DO 121 IJ=I,IP1

MNP(IJ)=MNP(IJ)-1

121 CONTINUE

105 CONTINUE

cwrite(6,66025) kis

66025format(1x,'faca posle 105 kis=',i3)

M=M-IP1

cwrite(*,*) m,kis

cpause 'm,kis'

104 CONTINUE

cwrite(*,*) kis

cpause '1'

GO TO 118

106 CONTINUE

cwrite(6,66034) kis

66034format(1x,'faca posle 106 kis=',i3)

M=M-1

118 CONTINUE

cwrite(*,*) kis

cpause '2'

cwrite(6,66013) kis

66013format(1x,' faca posle 118 kis=',i3)

c-формирование матрицы факторных нагрузок после уплотнения строк

c-элементы матрицы расположены по столбцам,верно!

DO 125 I=2,K1

MNE=(I-1)*M

MPE=(I-1)*MST+1

MKE=MPE+M-1

DO 125 IPE=MPE,MKE

MNE=MNE+1

V(MNE)=V(IPE)

125 CONTINUE

CALL VARMX(M,K1,V,NC,TV,H,F,D)

WRITE(6,9232)K1

cWRITE(6,502)

CALL GF(M,K1,V)

CALL GRPK(M,K1,V,MNPY,IP1,MNPI,KIS)

2998 CONTINUE

727 CONTINUE

5152 CONTINUE

cWRITE(6,7026)

RETURN

END

SUBROUTINE GF(M,K1,V)

DIMENSION MD(5)

dimension s1(151),s2(151)

DIMENSION V(1)

JK=K1*M

306 FORMAT(' ',5F10.3)

IN=1

IR=M

401 format(' ','фактор',i5)

402 FORMAT(' ','столбец матрицы нагрузок',I5)

DO 102 I=1,5

MD(I)=I

102 CONTINUE

105 FORMAT(' ',':',5('фaктор',I2,':'))

106 FORMAT(' ',51('-'))

107 FORMAT(' ',':',5(F9.4,':'))

c o

cWRITE(6,106)

cWRITE(6,105)MD

cWRITE(6,106)

DO 765 J=1,K1

S2(J)=0.

765 CONTINUE

DO 108 I=1,M

cWRITE(6,107)(V(J),J=I,JK,M)

L=0

DO 701 J=I,JK,M

L=L+1

S1(L)=V(J)**2

S2(L)=S2(L)+V(J)**2

701 CONTINUE

702 FORMAT(1X,'квадраты нагрузок')

703 FORMAT(1X,5F10.3)

cWRITE(6,106)

108 continue

WRITE(6,706)

706 FORMAT(/9X,'Таблица-Сумма квадратов нагрузок по факторам.'/)

write(6,1201)

write(6,1200)

write(6,1201)

1201format(10x,40('-'))

1200format(10x,'|НОМЕР ФАКТОРА|СУММА КВАДРАТОВ НАГРУЗОК|')

do 1704 j=1,k1

WRITE(6,1703) j,S2(J)

write(6,1201)

1704continue

1703format(10x,'|',5x,i3,5x,'|',9x,f6.3,9x,'|')

RETURN

END

subroutine corre (n,m,io,x,xbar,std,rx,r,b,d,t)

dimension x(1),xbar(1),std(1),rx(1),r(1),b(1),d(1),t(1)

do 100 j=1,m

b(j)=0.0

100 T(J)=0.0

k=(m*m+m)/2

do 102 i=1,k

102 R(I)=0.0

fn=n

l=0

if(io) 105,127,105

105 DO 108 J=1,M

do 107 i=1,n

l=l+1

107 T(J)=T(J)+X(L)

xbar(j)=t(j)

108 T(J)=T(J)/FN

do 115 i=1,n

jk=0

l=i-n

do 110 j=1,m

l=l+n

d(j)=x(l)-t(j)

110 B(J)=B(J)+D(J)

do 115 j=1,m

do 115 k=1,j

jk=jk+1

115 r(jk)=r(jk)+d(j)*d(k)

go to 205

127 IF(N-M) 130,130,135

130 KK=N

go to 137

135 KK=M

137 DO 140 I=1,KK

c call data (m,d)

do 140 j=1,m

t(j)=t(j)+d(j)

l=l+1

140 RX(L)=D(J)

fkk=kk

do 150 j=1,m

xbar(j)=t(j)

150 T(J)=T(J)/FKK

l=0

do 180 i=1,kk

jk=0

do 170 j=1,m

l=l+1

170 D(J)=RX(L)-T(J)

do 180 j=1,m

b(j)=b(j)+d(j)

do 180 k=1,j

jk=jk+1

180 R(JK)=R(JK)+D(J)*D(K)

if(n-kk) 205,205,185

185 KK=N-KK

do 200 i=1,kk

jk=0

c call data(m,d)

do 190 j=1,m

xbar(j)=xbar(j)+d(j)

d(j)=d(j)-t(j)

190 B(J)=B(J)+D(J)

do 200 j=1,m

do 200 k=1,j

jk=jk+1

200 R(JK)=R(JK)+D(J)*D(K)

205 JK=0

do 210 j=1,m

xbar(j)=xbar(j)/fn

do 210 k=1,j

jk=jk+1

210 R(JK)=R(JK)-B(J)*B(K)/FN

jk=0

do 220 j=1,m

jk=jk+j

220 STD(J)=SQRT(ABS(R(JK)))

do 230 j=1,m

do 230 k=j,m

jk=j+(k*k-k)/2

l=m*(j-1)+k

rx(l)=r(jk)

l=m*(k-1)+j

rx(l)=r(jk)

if(std(j)*std(k)) 225,222,225

222 R(JK)=0.0

go to 230

225 R(JK)=R(JK)/(STD(J)*STD(K))

230 CONTINUE

fn=sqrt(fn-1.0)

do 240 j=1,m

240 std(j)=std(j)/fn

l=-m

do 250 i=1,m

l=l+m+1

250 B(I)=RX(L)

return

end

SUBROUTINE EIGEN(A,R,N,MV) EIGE0010

DIMENSION A(1),R(1) EIGE0020

5 RANGE=1.0E-6 EIGE0030

IF(MV-1) 10,25,10 EIGE0040

10 IQ=-N EIGE0050

DO 20 J=1,N EIGE0060

IQ=IQ+N EIGE0070

DO 20 I=1,N EIGE0080

IJ=IQ+I EIGE0090

R(IJ)=0.0 EIGE0100

IF(I-J) 20,15,20 EIGE0110

15 R(IJ)=1.0 EIGE0120

20 CONTINUE EIGE0130

25 ANORM=0.0 EIGE0140

DO 35 I=1,N EIGE0150

DO 35 J=I,N EIGE0160

IF(I-J) 30,35,30 EIGE0170

30 IA=I+(J*J-J)/2 EIGE0180

ANORM=ANORM+A(IA)*A(IA) EIGE0190

35 CONTINUE EIGE0200

IF(ANORM) 165,165,40 EIGE0210

40 ANORM=1.414*SQRT(ANORM) EIGE0220

ANRMX=ANORM*RANGE/FLOAT(N) EIGE0230

IND=0 EIGE0240

THR=ANORM EIGE0250

45 THR=THR/FLOAT(N) EIGE0260

50 L=1 EIGE0270

55 M=L+1 EIGE0280

60 MQ=(M*M-M)/2 EIGE0290

LQ=(L*L-L)/2 EIGE0300

LM=L+MQ EIGE0310

62 IF( ABS(A(LM))-THR) 130,65,65 EIGE0320

65 IND=1 EIGE0330

LL=L+LQ EIGE0340

MM=M+MQ EIGE0350

X=0.5*(A(LL)-A(MM)) EIGE0360

68 Y=-A(LM)/ SQRT(A(LM)*A(LM)+X*X) EIGE0370

IF(X) 70,75,75 EIGE0380

70 Y=-Y EIGE0390

75 SINX=Y/ SQRT(2.0*(1.0+( SQRT(1.0-Y*Y)))) EIGE0400

SINX2=SINX*SINX EIGE0410

78 COSX= SQRT(1.0-SINX2) EIGE0420

COSX2=COSX*COSX EIGE0430

SINCS =SINX*COSX EIGE0440

ILQ=N*(L-1) EIGE0450

IMQ=N*(M-1) EIGE0460

DO 125 I=1,N EIGE0470

IQ=(I*I-I)/2 EIGE0480

IF(I-L) 80,115,80 EIGE0490

80 IF(I-M) 85,115,90 EIGE0500

85 IM=I+MQ EIGE0510

GO TO 95 EIGE0520

90 IM=M+IQ EIGE0530

95 IF(I-L) 100,105,105 EIGE0540

100 IL=I+LQ EIGE0550

GO TO 110 EIGE0560

105 IL=L+IQ EIGE0570

110 X=A(IL)*COSX-A(IM)*SINX EIGE0580

A(IM)=A(IL)*SINX+A(IM)*COSX EIGE0590

A(IL)=X EIGE0600

115 IF(MV-1) 120,125,120 EIGE0610

120 ILR=ILQ+I EIGE0620

IMR=IMQ+I EIGE0630

X=R(ILR)*COSX-R(IMR)*SINX EIGE0640

R(IMR)=R(ILR)*SINX+R(IMR)*COSX EIGE0650

R(ILR)=X EIGE0660

125 CONTINUE EIGE0670

X=2.0*A(LM)*SINCS EIGE0680

Y=A(LL)*COSX2+A(MM)*SINX2-X EIGE0690

X=A(LL)*SINX2+A(MM)*COSX2+X EIGE0700

A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) EIGE0710

A(LL)=Y EIGE0720

A(MM)=X EIGE0730

130 IF(M-N) 135,140,135 EIGE0740

135 M=M+1 EIGE0750

GO TO 60 EIGE0760

140 IF(L-(N-1)) 145,150,145 EIGE0770

145 L=L+1 EIGE0780

GO TO 55 EIGE0790

150 IF(IND-1) 160,155,160 EIGE0800

155 IND=0 EIGE0810

GO TO 50 EIGE0820

160 IF(THR-ANRMX) 165,165,45 EIGE0830

165 IQ=-N EIGE0840

DO 185 I=1,N EIGE0850

IQ=IQ+N EIGE0860

LL=I+(I*I-I)/2 EIGE0870

JQ=N*(I-2) EIGE0880

DO 185 J=I,N EIGE0890

JQ=JQ+N EIGE0900

MM=J+(J*J-J)/2 EIGE0910

IF(A(LL)-A(MM)) 170,185,185 EIGE0920170 X=A(LL) EIGE0930

A(LL)=A(MM) EIGE0940

A(MM)=X EIGE0950

IF(MV-1) 175,185,175 EIGE0960

175 DO 180 K=1,N EIGE0970

ILR=IQ+K EIGE0980

IMR=JQ+K EIGE0990

X=R(ILR) EIGE1000

R(ILR)=R(IMR) EIGE1010

180 R(IMR)=X EIGE1020

185 CONTINUE EIGE1030

RETURN EIGE1040

END EIGE1050

SUBROUTINE LOAD (M,K,R,V)

DIMENSION R(1),V(1)

L=0

JJ=0