Пи решении систем уравнений с матрицами коэффициентов, близкими к вырожденным, можно воспользоваться целочисленной арифметикой. Эта идея возникла у меня очень давно, гораздо раньше чем появились ПиСи и Си и Си++. На БЭСМ-6 был реализован алгоритм на Фортране. Листая свои записи я вспомнил о нем и проверил его на современном Inel Fortran под Visual Studio.NET-2008. Оказалось все работает. Думаю что это может пригодиться, тем более что сейчас многие вадеют программированием под MS Visual Studio.NET.
Вот текст:
Код:
S U B R O U T I N E LUFACT ( N , LU , IPS , SCALES , IER )
I M P L I C I T I N T E G E R * 4 ( A - Z )
C Подпрограмма реализует LU разложение исходной матрицы
C N - размерность системы уравнений
C LU- матрица коэффициентов
C IPS - вспомогательный массив
C SCALES - вспомогательный массив
C IER - индикатор ошибки
D I M E N S I O N IPS( N ) , SCALES( 2, N ) , ROWNRM( 2 )
* , BIG( 2 ) , ISIZE( 2 ), PIVOT( 2 ), EM( 2 ), PC( 2 )
* , LU( 2, N, N )
C Начальные установки и поиск ведущего элемента
IER = 0
D O 3 I = 1 , N
IPS( I ) = I
ROWNRM( 1 ) = 0
ROWNRM( 2 ) = 1
D O 1 J = 1 , N
IF ( ROWNRM( 1 )*LU( 2 , I , J ) .GE.
* ROWNRM( 2 )*IABS(LU( 1 , I , J ) ) ) GO TO 1
ROWNRM( 1 ) = IABS ( LU( 1 , I , J ) )
ROWNRM( 2 ) = LU( 2 , I , J )
1 C O N T I N U E
C Проверка ведущего элемента
IF ( ROWNRM( 1 ) ) 2 , 1001 , 2
2 SCALES( 1 , I ) = ROWNRM( 2 )
SCALES( 2 , I ) = ROWNRM( 1 )
3 C O N T I N U E
C Факторизация
NM1 = N-1
D O 9 K = 1 , NM1
BIG( 1 ) = 0
BIG( 2 ) = 1
DO 4 I = K , N
IP = IPS( I )
ISIZE( 1 ) = IABS ( LU( 1 , IP , K ) )*SCALES( 1 , IP )
ISIZE( 2 ) = LU( 2 , IP , K )*SCALES( 2 , IP )
C A L L SOKRAS ( ISIZE( 1 ) , ISIZE( 2 ) )
IF ( ISIZE( 1 )*BIG( 2 ) .LE. ISIZE( 2 )*BIG( 1 ) ) GO TO 4
BIG( 1 ) = ISIZE( 1 )
BIG( 2 ) = ISIZE( 2 )
PIVROW = I
4 C O N T I N U E
C Проверка на ошибку
IF ( BIG( 1 ) .EQ. 0 ) GO TO 1002
C Перестановки
IF ( PIVROW .EQ. K ) GO TO 5
J = IPS( K )
IPS( K ) = IPS( PIVROW )
IPS( PIVROW ) = J
5 KP = IPS( K )
PIVOT( 1 ) = LU( 2 , KP , K )
IF ( LU( 1 , KP , K ) ) 6 , 1002 , 7
6 PIVOT( 1 ) = -PIVOT( 1 )
7 PIVOT( 2 ) = IABS(LU( 1 , KP , K ) )
LU( 1 , KP , K) = PIVOT( 1 )
LU( 2 , KP , K) = PIVOT( 2 )
KP1 = K+1
D O 8 I = KP1 , N
IP = IPS( I )
EM( 1 ) = LU( 1 , IP , K )*PIVOT( 1 )
EM( 2 ) = LU( 2 , IP , K )*PIVOT( 2 )
C A L L SOKRAS ( EM( 1 ) , EM( 2 ) )
LU( 1 , IP , K ) = EM( 1 )
LU( 2 , IP , K ) = EM( 2 )
D O 8 J = KP1 , N
PC( 1 ) = EM( 1 )*LU( 1 , KP , J )
PC( 2 ) = EM( 2 )*LU( 2 , KP , J )
C A L L SOKRAS ( PC( 1 ) , PC( 2 ) )
LU( 1 , IP , J ) = LU( 1 , IP , J )*PC( 2 )-
* LU( 2 , IP , J )*PC( 1 )
LU( 2 , IP , J ) = LU( 2 , IP , J )*PC( 2 )
C A L L SOKRAS ( LU( 1 , IP , J ) , LU( 2 , IP , J ) )
8 CONTINUE
9 CONTINUE
C Перестановки с новым ведущим элементом
IP = IPS( N )
PIVOT( 1 ) = LU( 1 , IP , N )
PIVOT( 2 ) = LU( 2 , IP , N )
IF( PIVOT( 1 ) ) 10 , 1002 , 11
10 PIVOT( 2 ) = -PIVOT( 2 )
11 LU( 1 , IP , N ) = PIVOT( 2 )
LU( 2 , IP , N ) = IABS( PIVOT( 1 ) )
R E T U R N
C Блок регистрации ошибок
1001 IER = 1
R E T U R N
C142
1002 IER = 2
R E T U R N
E N D
S U B R O U T I N E LUSOLV ( N , LU , B , X , IPS )
I M P L I C I T I N T E G E R * 4 ( A - Z )
C Подпрограмма реализует решение системы ур-ний
C после LU разложения исходной матрицы
C N - размерность системы уравнений
C LU- матрица коэффициентов
C IPS - вспомогательный массив
C B - вектор правых частей
C Х - вектор корней системы
C
D I M E N S I O N B( 2 , N ) , X( 2 , N ) , IPS( N )
* , LU( 2 , N , N ) , ISUM( 2 ) , PC( 2 )
NP1 = N+1
C144
IP = IPS( 1 )
X( 1 , 1 ) = B( 1 , IP )
X( 2 , 1 ) = B( 2 , IP )
D O 2 I = 2 , N
IP = IPS( I )
IM1 = I-1
ISUM( 1 ) = 0
ISUM( 2 ) = 1
D O 1 J = 1 , IM1
PC( 1 ) = LU( 1 , IP , J )*X( 1 , J )
PC( 2 ) = LU( 2 , IP , J )*X( 2 , J )
C A L L SOKRAS ( PC( 1 ) , PC( 2 ) )
ISUM( 1 ) = ISUM( 1 )*PC( 2 )+ISUM( 2 )*PC( 1 )
ISUM( 2 ) = ISUM( 2 )*PC( 2 )
1 C A L L SOKRAS ( ISUM( 1 ) , ISUM( 2 ) )
X( 1 , I ) = B( 1 , IP )*ISUM( 2 )-B( 2 , IP )*ISUM( 1 )
X( 2 , I ) = B( 2 , IP )*ISUM( 2 )
2 C A L L SOKRAS ( X( 1 , I ) , X( 2 , I ) )
IP = IPS( N )
X( 1 , N ) = X( 1 , N )*LU( 1 , IP , N )
X( 2 , N ) = X( 2 , N )*LU( 2 , IP , N )
C A L L SOKRAS ( X( 1 , N ) , X( 2 , N ) )
C145
D O 4 IBACK = 2 , N
I = NP1-IBACK
IP = IPS( I )
IP1 = I+1
ISUM( 1 ) = 0
ISUM( 2 ) = 1
D O 3 J = IP1 , N
PC( 1 ) = LU( 1 , IP , J )*X( 1 , J )
PC( 2 ) = LU( 2 , IP , J )*X( 2 , J )
C A L L SOKRAS ( PC( 1 ) , PC( 2 ) )
ISUM( 1 ) = ISUM( 1 )*PC( 2 )+ISUM( 2 )*PC( 1 )
ISUM( 2 ) = ISUM( 2 )*PC( 2 )
3 C A L L SOKRAS( ISUM( 1 ) , ISUM( 2 ) )
PC( 1 ) = X( 1 , I )*ISUM( 2 )-X( 2 , I )*ISUM( 1 )
PC( 2 ) = X( 2 , I )*ISUM( 2 )
C A L L SOKRAS ( PC( 1 ) , PC( 2 ) )
X( 1 , I ) = PC( 1 )*LU( 1 , IP , I )
X( 2 , I ) = PC( 2 )*LU( 2 , IP , I )
4 C A L L SOKRAS( X( 1 , I ) , X( 2 , I ) )
R E T U R N
E N D
S U B R O U T I N E SOKRAS ( IA , IB )
I M P L I C I T I N T E G E R * 4 ( A - Z )
C Подпрограмма сокращения дробей
C IA - числитель (он же носитель знака числа)
C IB - знаменатель
C
IF (IA.EQ.0) IB=1
IF ( IB .LE. 0 ) S T O P 9999
IK = MIN ( IABS( IA ) , IABS ( IB ) )
I F ( IK .LE. 1 ) R E T U R N
99 C O N T I N U E
LCH = IA/IK
IF ( LCH*IK .NE. IA ) GO TO 1
LZN = IB/IK
IF ( LZN*IK .NE. IB ) GO TO 1
IA = LCH
IB = LZN
GO TO 99
1 C O N T I N U E
IK = IK-1
IF ( IK .GT. 1 ) GO TO 99
R E T U R N
E N D
