Eugeen
Member | Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору В старое (доброе время) когда не было четверной точности, а лучшей машиной была БЭСМ-6, я написал небольшой код для точных решений систем линейных уравнений (в целочисленной арифметике на обыкновенных дробях). Может кому пригодится: 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 ) C135 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 C136 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 C137 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 C138 IF ( BIG( 1 ) .EQ. 0 ) GO TO 1002 C13 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 C140 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 C141 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 I M P L I C I T I N T E G E R * 4 ( A - Z ) C C Тестовая программа C PARAMETER(N=3) I N T E G E R * 4 LU( 2,N,N) DIMENSION IPS( N ) , SCALES(2,N) C143 DIMENSION B(2,N),X(2,N) CALL BACKW1(N,LU,B,SCALES,IPS,X) STOP 10 END SUBROUTINE BACKW1(N,LU,B,SCALES,IPS,X,IERR) I M P L I C I T I N T E G E R * 4 ( A - H , O - Z ) D I M E N S I O N IPS( N ) , SCALES( N , 2 ) * , B( N , 2 ) , X ( N , 2 ) , LU( N , N , 2) DO 1 I=1,N B(I,1)=0 B(I,2)=1 DO 1 J=1,N LU(J,I,1)=1 1 LU(J,I,2)=1 B(2,1)=1 DO 2 I=1,N DO 2 J=2,N 2 LU(I,J,1)=(1-I)**(J-1) CALL DECOMP (N,N,LU,IPS,SCALES,IERR) CALL SOLVE(N,N,LU,B,X,IPS) WRITE (*,10) (X(I,1),X(I,2), I=1,N) RETURN 10 FORMAT(1X,1P,I5,'/',I5) END |