C C Программа PseudsysccmSolver - решения систем линейных уравнений AX=Y C со ``слабо'' вырожденной (у A имеется лишь одно нулевое собственное C значение) матрицей A (метод критических компонент [88]) C C Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN C Обpащение: CALL PseudsysccmSolver(M,A,N,INF,IM,B,R,C,P) C Входные данные: Блок F499_RCONST содержит вещественные константы C e_8,e_0,e_1,e_2=e_1/beta определенные в INIT_CONST C (или заранее известные) C M,N - (integer) размерности прямоугольной матpицы A; C INF - (integer): C INF=1, если pешается система с верхнедвухдиагональной C матpицей; C INF=-1, если pешается система с нижнедвухдиагональной C матpицей; C INF=2, если pешается (без масштабиpования) система с C заполненной A матpицей; C INF=-2, если pешается (с масштабиpованием) система с C заполненной A матpицей; C A - (real*8) двумеpный массив (размеpность массива A: C [M,N], если на входе |INF|=2; [1,1] (или A - C вещественная пеpеменная), если на входе |INF|=1), C котоpый содеpжит исходную матpицу A. Пpи этом матpицы C A запоминаются в виде (****); C IM - (integer) число правых частей; C B - (real*8) двумеpный массив. Имеет размерность C [max(M,N),IM]. При этом в B хранится матрица правых C частей в виде (****); C R - (real*8) одномеpный pабочий массив pазмеpности 2М. C В нем на входе содеpжится (пpи INF=-1 и INF=1) C двухдиагональная матpица в виде (*); C C,P - (real*8) pабочые массивы pазмеpности 2М и 3М C соответственно. C Выходные данные: INF - (integer): C INF=0 - ноpмальное завеpшение pаботы подпpогpамм; C INF=-1 - исходная матpица выpождена; C A - (real*8) массив содеpжит матpицы P и Q [10] в C в разложении PAQ=C_2; C B - (real*8) массив содеpжит матpицу pешений в виде (****) C SUBROUTINE PseudsysccmSolver(M,A,N,INF,IM,B,R,C,P) IMPLICIT REAL*8 (A - H, O - Z) DIMENSION A(M,*),B(*),R(*),C(*),P(*) INTEGER EAB,EX,MSCL COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN C call Init_Const ! Sap: to hide this actin from User IF(M.LE.0.OR.N.LE.0) write(*,100)M,N M0=INF IM0=IABS(M0) MXM=MAX0(M,N) MNM=MIN0(M,N) M1=MNM-1 IF(MXM.LE.0) RETURN INF=-1 IF(IM0.EQ.1) GOTO 3 * IF(M0.GT.0) GOTO 1 IF(MSCL(0,M,N,A,EAB).EQ.0) RETURN 1 CALL BTD(IM0,M,N,A,R) EX=0 IF(M0.GT.0) GOTO 2 IF(MSCL(0,2,M,R,EX).EQ.0) RETURN EAB=EAB+EX IF(MSCL(0,MXM,IM,B,EX).EQ.0) CONTINUE EX=EX-EAB 2 CALL BTDBK(IM0,M,N,A,1,1,IM,B) 3 IF(DABS(R(MNM*2)).LT.UNDFLO) GOTO 4 * IF ((M0.EQ.1).OR.((M0.EQ.-1).AND.(M.GE.N))) INF=1 DET=0D0 CALL Lin2dsysccmSolver(MNM,R,M,INF,I,JM0,B,DET) IF(IM0.EQ.1) GOTO 18 GOTO 13 4 IF(M0.EQ.-1) GOTO 8 DO 5 I=1,MNM-2 C(I*2)=R(I*2)*R(I*2)+R(I*2+1)*R(I*2+1) C(I*2+1)=R(I*2+2)*R(I*2+1) 5 CONTINUE C(M1*2)=R(M1*2)*R(M1*2)+R(MNM*2-1)*R(MNM*2-1) INF=1 DET=0D0 CALL Lin3dsysccmSolver(M1,C,M,INF,IM,0,B,P,DET) DO 7 L=1,IM B((L-1)*M+MNM)=R(MNM*2-1)*B((L-1)*M+M1) DO 6 I=M1,2,-1 B((L-1)*M+I)=R(I*2-1)*B((L-1)*M+I-1)+R(I*2)*B((L-1)*M+I) 6 CONTINUE B((L-1)*M+1)=R(2)*B((L-1)*M+1) 7 CONTINUE GOTO 13 C 8 DO 9 I=MNM,3,-1 C(I*2-2)=R(I*2)*R(I*2)+R(I*2-1)*R(I*2-1) C(I*2-3)=R(I*2-2)*R(I*2-1) 9 CONTINUE C(2)=R(4)*R(4)+R(3)*R(3) DO 10 J=1,IM DO 10 I=1,M1 10 B((J-1)*M+I)=B((J-1)*M+I+1) INF=1 DET=0D0 CALL Lin3dsysccmSolver(M1,C,M,INF,IM,0,B,P,DET) DO 12 L=1,IM B((L-1)*M+1)=R(3)*B((L-1)*M+1) DO 11 I=2,M1 B((L-1)*M+I)=R(I*2+1)*B((L-1)*M+I)+R(I*2)*B((L-1)*M+I-1) 11 CONTINUE B((L-1)*M+MNM)=R(2*MNM)*B((L-1)*M+MNM) 12 CONTINUE C 13 IF(M.EQ.N) GOTO 17 NI=N*(I-1) IF(M.GT.N) THEN L=1 DO 14 K=N,NI,N L=L+1 DO 14 J=1,N 14 B(K+J)=B((L-1)*M+J) ELSE L=I+1 DO 15 K=NI,N,-N L=L-1 DO 15 J=M,1,-1 B(K+J)=B((L-1)*M+J) 15 B((L-1)*M+J)=0D0 DO 16 J=M+1,N 16 B(NI+J)=0D0 ENDIF 17 CALL BTDBK(4,M,N,A,-1,1,IM,B) IF(M0.GT.0) GOTO 18 IF(MSCL(1,MXM,IM,B,EX).EQ.0) CONTINUE * 18 INF=0 RETURN 100 FORMAT(' M OR N LESS THEN 0 ON PROGRAM PseudsysccmSolver: M=', * I2,' N=',I2) END C C Программы полученные из библиотеки LINA [3] (А.Н.Малышев, 1990) C Программа BTD C C Приведение вещественной матрицы к двухдиагональному виду или C симметричной вещественной матрицы к трехдиагональному виду C преобразованиями отражения Хаусхолдера C C Обpащение: CALL BTD(I23,M,N,A,B) C Входные данные: C I23 - (integer) равно 2, если приведение к C двухдиагональному виду; равно 3, если приведение C к трехдиагональному виду; C A - (real*8) иходная вещественная матрица; C M - (integer) число строк матрицы A; C N - (integer) число столбцов матрицы A (при I23=3 C должно быть N=M). C Выходные данные: C B - (real*8) массив длины 2*L, где L=min(M,N). C B(2),B(4),...,B(2*L) - элементы главной диагонали; C B(3),B(5),...,B(2*L-1) - элементы побочной диагонали; C A - (real*8) в этом массиве упакованы векторы нормалей C отражений; C SUBROUTINE BTD(I23,M,N,A,B) INTEGER I23,M,N,INC,I,II,I1,L,LL,L1,J,K,K1,K2 REAL*8 A(*),B(*),RFLD,ZERO LOGICAL I2 DATA ZERO/0D0/ I2=I23.EQ.2 B(1)=ZERO K2=1 K=I23-1 L=M+1-K IF(L.LE.0) GOTO 5 LL=N I=1 II=M I1=1 J=1 INC=K IF(M.GE.N) GOTO 3 1 I1=II II=I I=I1 IF(I2) GOTO 2 IF(I.EQ.1) GOTO 2 B(J-1)=A(K1-1) I1=1 GOTO 4 2 L1=LL LL=L L=L1 3 J=J+INC B(J)=RFLD(L,A(K),I) K1=K K2=K+II L1=LL-1 4 LL=LL-1 IF(LL.LE.0) GOTO 5 K=K+II IF(L.GT.1) CALL REFL(L,A(K1),I1,L1,A(K2),I,II) GOTO 1 5 IF(.NOT.I2) B(N+N)=A(K2) RETURN END C C Программа BTDBK C C Применение преобразований отражений Хаусхолдера, вычисленных C в подпрограмме BTD, к вещественной матрице C C Обpащение: CALL BTDBK(I234,M,N,A,I,LR12,K,X) C Входные данные: C A - (real*8) матрица, содержащая векторов нормалей C отражений, упакованные подпрограммой BTD; C M - (integer) число строк матрицы A; C N - (integer) число столбцов матрицы A (при I23=3 C должно быть N=M). C I234 - (integer) равно 2, если применяются левые отражения C из процесса двухдиагонализации; равно 4, если C применяются правые отражения из процесса C двухдиагонализации; равно 3, если применяются левые C отражения из процесса трехдиагонализации; C I - (integer) равно 1, если отражения применяются в прямом C порядке; равно -1, если отражения применяются в C обратном порядке; C X - (real*8) исходная матрица, к которой применяются C преобразования отражения, упакованные в матрице A; C LR12 - (integer) равно 1, если отражения применяются к C к матрице X слева; равно 2, если отражения C применяются к матрице X справа; C K - (integer) при LR12=1 число столбцов, а при LR12=2 C число строк матрицы X. C Выходные данные: C X - (real*8) содержит результат применения преобразований C отражений к исходной матрице X. C SUBROUTINE BTDBK(I234,M,N,A,I,LR12,K,X) INTEGER I234,M,N,I,LR12,K REAL*8 A(*),X(*) INTEGER IA,IX,JX,L,LL,LLL,J,JJ,MP1I,IXI LOGICAL I4 I4=I234.EQ.4 IA=1 IF(I4) IA=M IX=K JX=1 IF(LR12.EQ.2) GO TO 1 IX=1 JX=M IF(I4) JX=N 1 L=MAX0(M,N) LL=MIN0(M,N) J=1 JJ=1 IF((I234.EQ.2).AND.(M.GE.N).OR.(I4.AND.(M.LT.N))) GOTO 2 L=LL-1 LL=L-1 J=J+IX JJ=JJ+IA GOTO 3 2 IF(M.EQ.N) LL=LL-1 3 IF(LL.LE.0) GOTO 6 MP1I=(M+1)*I IXI=IX*I IF(I.EQ.1) GOTO 4 LLL=LL-1 L=L-LLL J=J-IXI*LLL JJ=JJ-MP1I*LLL 4 DO 5 LLL=1,LL CALL REFL(L,A(JJ),IA,K,X(J),IX,JX) L=L-I J=J+IXI 5 JJ=JJ+MP1I 6 RETURN END C C Программа REFL C C Применение одного преобразования отражения Хаусхолдера к C вещественной матрице C C Обpащение: CALL REFL(L,P,IP,N,X,IX,JX) C Входные данные: C L - (integer) число компонент в нормали, определяющей C преобразование отражения; C P - (real*8) массив, содержащий компоненты P(1),P(1+IP), C ...,P(1+(L-1)*IP) вектора нормали отражения; C IP - (integer) интервал между компонентами вектора C нормали отражения в массиве P; C N - (integer) число столбцов в преобразуемой матрице; C X - (real*8) массив, содержащий элементы C C X(1), X(1+JX), ...,X(1+(N-1)*JX), C X(1+IX), X(1+IX+JX), ...,X(1+IX+(N-1)*JX), C : : : : C X(1+(L-1)*IX),X(1+(L-1)*IX+JX),...,X(1+(L-1)*IX+(N-1)*JX); C C IX - (integer) интервал внутри столбцов между компонентами C преобразуемой матрицы в массиве X; C JX - (integer) интервал внутри строк между компонентами C преобразуемой матрицы в массиве X. C Выходные данные: C X - (real*8) в элементах X(1+I*IX+J*JX), где I=0,1,...,L-1, C J=0,1,...,N-1, содержится преобразованная матрица C SUBROUTINE REFL(L,P,IP,N,X,IX,JX) INTEGER I,J,K,IPL,JXN,L,IP,N,IX,JX REAL*8 INPR,S,P(*),X(*) IPL=IP*L JXN=JX*N DO 1 J=1,JXN,JX S=INPR(L,P,IP,X(J),IX) K=J DO 1 I=1,IPL,IP X(K)=X(K)-S*P(I) 1 K=K+IX RETURN END C C Функция RFLD C C Построение преобразования отражения Хаусхолдера, которое заданный C вещественной вектор переводит в вектор, имеющий нулевые компоненты C со второй по последнюю включительно C C Обpащение: RFLD(L,P,IP) C Входные данные: C L - (integer) размерность заданного вектора; C P - (real*8) массив, содержащий компоненты P(1),P(1+IP), C ...,P(1+(L-1)*IP) заданного вектора, по которому C строится преобразование отражения; C IP - (integer) интервал между компонентами заданного C вектора в массиве P. C Выходные данные: C P - (real*8) в элементах P(1),P(1+IP),...,P(1+(L-1)*IP) C получены компоненты вектора нормали, определяющий C преобразования отражения; C RFLD - (real*8) первая компонента вектора, полученного при C отражении заданного вектора C REAL*8 FUNCTION RFLD(L,P,IP) INTEGER I,IPL,L,IP REAL*8 PM,T,ZERO,INPR,P(*) DATA ZERO/0D0/ IF(L.EQ.1) GO TO 5 IPL=IP*L PM=ZERO DO 1 I=1,IPL,IP 1 PM=DMAX1(PM,DABS(P(I))) IF(PM.LE.ZERO) GO TO 5 DO 2 I=1,IPL,IP 2 P(I)=P(I)/PM T=DSQRT(INPR(L,P,IP,P,IP)) IF(P(1).LE.ZERO) T=-T RFLD=-T*PM P(1)=P(1)+T PM=DSQRT(P(1)*T) DO 3 I=1,IPL,IP 3 P(I)=P(I)/PM 4 RETURN 5 RFLD=P(1) P(1)=ZERO GO TO 4 END C C Функция MSCL C C Вычисления масштабного разложения вещественной матрицы или C умножение вещественной матрицы на степень основания счисления C машинных вещественных чисел C C Обpащение: MSCL(I01,M,N,A,E) C Входные данные: C I01 - (integer) равно 0, если вычисляется масштабное C разложение; равно 1, если выполняется умножение C на степень основания счисления; C A - (real*8) исходная обрабатываемая матрица; C M - (integer) число строк матрицы A; C N - (integer) число столбцов матрицы A; C E - (integer) при I01=1 показатель степени основания C счисления. C Выходные данные: C A - (real*8) при I01=0 результат масштабирования исходной C матрицы A; при I01=1 равно произведению исходной C матрицы A на степень основания счисления с C показателем E; C E - (integer) при I01=0 экспонента в масштабном разложении C исходной матрицы A; C MSCL - (integer) равно -1, если умножение на степень C основания невозможно из-за переполнения (в этом C случае матрица A на выходе равна исходной); C равно 0, если на входе A=0; C равно 1 в остальных случаях. C INTEGER FUNCTION MSCL(I01,M,N,A,E) INTEGER MN,I,EE,E,I01,M,N REAL*8 AM,ZERO,GTEXP,CHEXP,A(*) INTEGER RADIX,MAXEXP,MINEXP,MANTSZ COMMON /F499_ICONST/RADIX,MAXEXP,MINEXP,MANTSZ DATA ZERO/0D0/ MN=M*N AM=ZERO DO 1 I=1,MN 1 AM=DMAX1(AM,DABS(A(I))) IF(AM.LE.ZERO) GO TO 7 AM=GTEXP(AM,EE) IF(I01.NE.0) GO TO 2 E=EE EE=-EE GO TO 3 2 IF(E.GT.MAXEXP-EE) GO TO 8 EE=E 3 IF(EE.EQ.0) GO TO 5 DO 4 I=1,MN 4 A(I)=CHEXP(A(I),EE) 5 MSCL=1 6 RETURN 7 IF(I01.EQ.0) E=0 MSCL=0 GO TO 6 8 MSCL=-1 GO TO 6 END C C Функция INPR C C Вычисления скалярного произведения вещественных векторов C C Обpащение: INPR(L,X,IX,Y,IY) C Входные данные: C L - (integer) число компонент у перемножаемых векторов; C X - (real*8) массив, содержащий в элементах X(1),X(1+IX), C ...,X(1+(L-1)*IX) компоненты первого вектора; C IX - (integer) интервал между компонентами первого C вектора в массиве X; C Y - (real*8) массив, содержащий в элементах Y(1),Y(1+IY), C ...,Y(1+(L-1)*IY) компоненты второго вектора; C IY - (integer) интервал между компонентами второго C вектора в массиве Y. C Выходные данные: C INPR - (real*8) равно X(1)*Y(1)+X(1+IX)*Y(1+IY)+... C +X(1+(L-1)*IX)+Y(1+(L-1)*IY) C REAL*8 FUNCTION INPR(L,X,IX,Y,IY) INTEGER I,J,IXL,L,IX,IY REAL*8 ZERO,X(*),Y(*) DATA ZERO/0D0/ INPR=ZERO IXL=IX*L J=1 DO 1 I=1,IXL,IX INPR=INPR+X(I)*Y(J) 1 J=J+IY RETURN END