c  ARYTMETYKA fl~ OBLICZEN IL. SKALARNYCH
c
c       n - rzad macierzy układu  rownan - ilosc rownan
c       nb = n*n
c       c(nb) - macierz układu  a(n,n) - po wykonaniu obliczen jednostkowa
c       d(np) - prawe strony ( m wektorow o dlugosci n) d(n,m) na wejsciu,
c               zapamietane wiersz po wierszu
c       m - ilosc prawych stron czyli ilosc układow rownan
c       np = n*m
c       d(np) - na wyjsciu ROZWIAZANIA d(n,m) zapamietane wiersz po wierszu
c       ierr - wskaznik bledu
c               ierr=0 nie bylo bledu
c               ierr=-1 macierz osobliwa

_________________________________________________________________________

        subroutine lgaus4(n,m,a,c,b,d,ierr)
c       NIE NISZCZY MACIERZY !!!
c  rozwiązanie m UKŁADOW ROWNAN LINIOWYCH METODA GAUSS'A Z KOLUMNOWYM
c     WYBOREM ELEMENTU PODSTAWOWEGO
c           a(n,n)*d(n,m)=b(n,m)        gdzie d(n,m)-niewiadome
c  ARYTMETYKA fl~ OBLICZEN IL. SKALARNYCH
c
c       n - rzad macierzy układu  rownan - ilosc rownan
c       nb = n*n
c       a(nb) - macierz układu  a(n,n) - po wykonaniu obliczen NIEZMIENONA
c       b(np) - prawe strony ( m wektorow o dlugosci n) b(n,m), zapamietane
c               wiersz po wierszu
c       c(nb) - robocza tablica dla a(nb)
c       m - ilosc prawych stron czyli ilosc układow rownan
c       np = n*m
c       d(np) - robocza tablica dla b(np) - na wyjsciu ROZWIAZANIA
c       ierr - wskaznik bledu
c               ierr=0 nie bylo bledu
c               ierr=-1 macierz osobliwa
__________________________________________________________________________
        subroutine aminusb(n,a,b,c)
c odejmowanie 2 macierzy o tych samych wymiarach
c
c wszystkie macierze w postaci wektorow
c c=c(n) , a=a(n) , b=b(n)
c zanotowane wiersz po wierszu
_________________________________________________________________________
        subroutine arazyb(n,m,l,a,b,c)
c OBLICZENIA ILOCZYNOW SKALARNYCH W ARYTMETYCE fl~
c mnozenie 2 macierzy prostokatnych
c c(n,l) = a(n,m) * b(m,l)
c
c wszystkie macierze w postaci wektorow
c c(n,l)=c(nc) , a(n,m)=a(na) , b(m,l)=b(nb)
c       nc=n*l ,   na=n*m     ,   nb=m*l
c zanotowane wiersz po wierszu
c
__________________________________________________________________________
        subroutine aplusb(n,a,b,c)
c dodawanie 2 macierzy o tych samych wymiarach
c
c wszystkie macierze w postaci wektorow
c c=c(n) , a=a(n) , b=b(n)
c zanotowane wiersz po wierszu
c
_________________________________________________________________________
        subroutine inversa1(n,c,b,ierr)
c               PROCEDURA NISZCZACA MACIERZ !!!
c               metoda Jordana
c       OBLICZANIE MACIERZY ODWROTNEJ DO MACIERZY c
c               PODWOJONA PRECYZJA real*8
c
c       n - rzad macierzy kwadratowej a
c       nb = n*n
c       c(nb) - macierz dana - na wyjsciu jednostkowa
c       b(nb) - na wyjsciu macierz odwrotna do c
c       ierr - wskaznik osobliwosci
c               ierr=0 macierz NIEosobliwa
c               ierr=-1 macierz OSObliwa
__________________________________________________________________________

        subroutine inversa2(n,a,c,b,ierr)
c               PROCEDURA NIENISZCZACA MACIERZY
c       OBLICZANIE MACIERZY ODWROTNEJ DO MACIERZY a
c               PODWOJONA PRECYZJA real*8
c
c       n - rzad macierzy kwadratowej a
c       nb = n*n
c       a(nb) - macierz dana - na wyjsciu nie zaburzona
c       b(nb) - na wyjsciu macierz odwrotna do a
c       c(nb) - robocza tablica dla a(nb)
c       ierr - wskaznik osobliwosci
c               ierr=0 macierz NIEosobliwa
c               ierr=-1 macierz OSObliwa
_________________________________________________________________________
        subroutine normae(na,a,euk)
        real*8 a(*),euk,u,p,v,s,z
c  NORMA EUKLIDESOWA MACIERZY KWADRATOWEJ a
c  obliczenia w arytmetyce fl~ (algorytm Gilla-Mollera)
c       n - rzad macierzy a
c       na=n*n
c       a(na) - macierz kwadratowa a
c       euk - norma euklidesowa macierzy a
c Dla macierzy jednostkowej wynosi sqrt(n)
__________________________________________________________________________
        subroutine okresla(n,a,c,e,d,ir,det,ierr)
c   BADANIE CZY MACIERZ JEST DODATNIO, SLABO DODATNIO LUB UJEMNIE
c                       OKRESLONA.
c
c       n - wymiar macierzy kwadratowej badanej a(n,n)
c       a(n*n) - macierz badana zapisana wiersz po wierszu
c       c(n*n) - macierz pomocnicza (po wykonaniu podprogramu zawiera
c                elementy macierzy L i U rozkladu na  czynniki trojkatne
c                o ile 'det' nie jest 0 )
c       e(n*n) - macierz pomocnicza
c       d(3*n) - na wejsciu moze zawierac prawe strony układu  rownan
c                d(1),...,d(n); o ile det nie rowny 0 to d(2n+1)...d(3n)
 

c                zawieraja elementy przekatnej glownej mac U z rozkladu
c                trojkatnego (uzyte do obliczenia wyznacznika)
c       ir(n) - tablica calkowita pomocnicza
c              Jesli ir(1)=-1 to nie jest obliczana wartosc wyznacznika 'det'
c               tylko podawany w 'det' jest jego znak
c       det - wyznacznik lub znak wyznacznika (patrz wyzej)
c       ierr - wskaznik wyniku dzialania podprogramu
c               ierr<0 macierz ujemnie okreslona
c               ierr=0 macierz nieujemnie okreslona lub osobliwa
c                       wyznacznik det=0 oznacza ze jest osobliwa
c               ierr>0 macierz dodatnio okreslona A > 0
_________________________________________________________________________
      subroutine thomas(n,a,b,c,d,x,bm,qm)
c rozwiązanie UKŁADU  ROWNAN Z MAC TROJPRZEKATNIOWA ALG. THOMAS'A
c
c       b1 c1                   = d1
c       a2 b2 c2                = d2
c          a3 b3 c3             = d3
c             a4 b4 c4          = d4
c
c       n - ilosc rownan
c
__________________________________________________________________________

        subroutine uwarun(n,a,c,b,uw,deti,ierr)
c       NIE NISZCZY MACIERZY !!!!
c  TESTOWANIE UWARUNKOWANIA MACIERZY KWADRATOWEJ a
c
c       n - rzad macierzy kwadratowej a(n,n)
c       na=n*n
c       a(na) - macierz a(n,n) - na wyjsciu niezaburzona
c       c(na) - macierz robocza
c       b(na) - macierz odwrotna do a(n,n)
c       uw - wskaznik uwarunkowania macierzy NORMAE(a)*NORMAE(b)
c       deti - norma iloczynu macierzy a i odwrotnej
c              powinna byc = sqrt(n)
c       ierr - wskaznik bledow :
c          = 0 ; macierz dobrze uwarunkowana
c          = 1 ; macierz zle lub potencjalnie zle uwarunkowana
c          =-1 ; macierz osobliwa-wykryte w algorytmie Jordana
c
----------------------------------------------------------------------------
        subroutine banach0(n,a,d,ierr)
c  ROZKLAD NA CZYNNIKI TROJKATNE MACIERZY a(n,n) SYMETRYCZNEJ I DODATNIO
c  OKRESLONEJ (NIEKONIECZNIE!) METODA Choleskiego-Banachiewicza
c    A(n,n) = L(n,n)*D(n,n)*L'(n,n)
c WSZYSTKIE TABLICE PAMIETANE WIERZS PO WIERZSU !!!
c       a(n,n)-dana macierz ( musi byc wypelniona tylko gorna polowa
c               wraz z diagonala - bo symetryczna),
c U W A G A :   jesli okaze sie ze
c               NIE jest SYMETRYCZNA (sprawdza tylko 2 ostanie wiersze
c               i kolumny) przepisze czesc gorna trojkatna do dolnej i
c               w ten sposob uczyni ja trojkatna !!!!
c       n - rzad macierzy a(n,n)
c  UWAGA:
c       po wykonaniu macierz L(n,n) ( tylko jej czesc dolna bez diagonali
c       bo diagonala to jedynki) zapisywana jest do macierzy a pod przekatna
c       tzn : L(i,j)=a(i,j) i=1,n;j=1,i-1
c       Diagonala macierzy L(n,n) zapisana do d(n) to znaczy ze:
c       tzn : L(i,i)=d(i)  i=1,n
c    ierr - wskaznik bledow
c       ierr < 0  - w trakcie rozkladu napotkano zero w d(i)
c       ierr > 0  - nie bylo bledow
c       ierr = 0  - macierz a(n,n) jest niesymetryczna ale zostala
c                   zsymetryzowana (tzn gorny trojkat zostal przepisany
c                   w dolny)
_________________________________________________________________________
        subroutine banach1(n,a,d,ierr)
c  rozwiązanie UKŁADU  ROWNAN Z MACIERZA SYMETRYCZNA I DODATNIO OKRESLONA
c   METODA ROZKLADU NA CZUNNIKI TROJKATNE Choleskiego-Banachiewicza
c               A(n,n) = L(n,n)*D(n,n)*L'(n,n)
c Niszczy prawe strony i dolna polowke macierzy A(n,n) (bez przekatnej)
c w dolnej polowce jest macierz L(n,n) bez diagonali do dalszego
c  wykorzystania
c a diagonala w tablicy d poczawszy od d(n+1)...d(2*n)
c
c       a(n*n) - macierz układu  symetryczna i A>0 , wystarczy tylko gorna
c               polowa
c       d(2*n) - tablica zawierajaca : d(1)...d(n) -prawe strony na WEJSCIU
c                              d(1)...d(n) - rozwiązanie na WYJSCIU
c                            d(n+1)...d(2*n)-diagonale D(n,n) na WYJSCIU
c       ierr - wskaznik bledu
c               ierr < 0 wystapil blad przy rozkladzie na czynniki trojkatne
 

c               ierr > 0 nie bylo bledu
        subroutine banach2(n,a,d)
c  rozwiązanie UKŁADU  ROWNAN Z MACIERZA SYMETRYCZNA I DODATNIO OKRESLONA
c   METODA ROZKLADU NA CZYNNIKI TROJKATNE Choleskiego-Banachiewicza
c               A(n,n) = L(n,n)*D(n,n)*L'(n,n)
c
c K O R Z Y S T A   Z WCZESNIEJ DOKONANEGO ROZKLADU PROGRAMEM banach1()
c
c       a(n*n) - macierz układu  symetryczna i A>0 :
c                       w gornej polowce macierz układu
c               w dolnej polowce macierz trojkatna dolna L z rozkladu
c                               wykonanego programem  banach1()
c       d(2*n) - tablica zawierajaca : d(1)...d(n) -prawe strony na WEJSCIU
c                                      d(1)...d(n) - rozwiązanie na WYJSCIU
c                       d(n+1)...d(2*n)-diagonale D(n,n) z poprzednio
c                                    wykonanego programu banach1()
__________________________________________________________________________
        subroutine banach3(n,m,a,d,ierr)
c  rozwiązanie m UKŁADOW ROWNAN Z TA SAMA MACIERZA SYMETRYCZNA I DODATNIO
c                               OKRESLONA
c   METODA ROZKLADU NA CZUNNIKI TROJKATNE Choleskiego-Banachiewicza
c               A(n,n) = L(n,n)*D(n,n)*L'(n,n)
c Niszczy prawe strony i dolna polowke macierzy A(n,n) (bez przekatnej)
c w dolnej polowce jest macierz L(n,n) bez diagonali do dalszego
c wykorzystania
c a diagonala w tablicy d poczawszy od d(m*n+1)...d(m*n+n)
c
c       n - rzad macierzy a(n,n)
c       m - ilosc prawych stron ( ilosc układow rownoczesnie rozwiazywanych)
c       a(n*n) - macierz układu  symetryczna i A>0 , wystarczy tylko gorna
c               polowa
c       d(m*n+n) - tablica : d(1)...d(n) -prawe strony na WEJSCIU 1 układu
c                  d(n+1)...d(2*n)-prawe strony na WEJSCIU 2 układu
c                    .          .
c             d((m-1)*n+1)...d(n*m)-prawe strony na WEJSCIU m-tego układu
 

c                       ----------------------------------------------------
c                d(1)...d(n) - rozwiązanie na WYJSCIU 1 układu
c              d(n+1)...d(2*n)-prawe strony na WYJSCIU 2 układu
c                    .          .
c             d((m-1)*n+1)...d(n*m)-prawe strony na WYJSCIU m-tego układu
 

c                       -------------------------------------------------
c                   d(n*m+1)...d(n*m+n)-diagonale D(n,n) na WYJSCIU
c       ierr - wskaznik bledu
c               ierr < 0 wystapil blad przy rozkladzie na czynniki trojkatne
 

c               ierr > 0 nie bylo bledu
__________________________________________________________________________
        subroutine arazybt(n,l,m,a,b,c)
c Oblicza iloczyn macierzy :
c               a * transpose{b}
c  gdzie macierze :
c       a(n,l) - pamietana wiersz po wierszu
c       b(m,l) - pamietana wiersz po wierszu ale nie transponowana !
c       c(n,m) - pamietana wiersz po wierszu
        real*8 a(*),b(*),c(*)
         subroutine atrazyb(n,l,m,a,b,c)
c Oblicza iloczyn macierzy :
c               transpose{a} * b
c  gdzie macierze :
c       a(l,n) - pamietana wiersz po wierszu ale nie transponowana
c       b(l,m) - pamietana wiersz po wierszu
c       c(n,m) - pamietana wiersz po wierszu
        real*8 a(*),b(*),c(*)
__________________________________________________________________________
        subroutine skrazya(sk,nn,a)
c Oblicza iloczyn macierzy :
c               sk * a
c  gdzie macierze :
c       a(nn) - pamietana wiersz po wierszu dowolna macierz
c       sk - sklar
        real*8 a(*),sk
__________________________________________________________________________
        subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1)
        integer qx,q1,q2,q3,q4,q5,q6,q7,q8,q9
c
        integer i,j,k,l,m,ii,i1,kk,k1,ll,l1,mn,nm,its,ierr
c       real*8 a(nm,n),w(n),u(nm,n),v(nm,n),rv1(n)
        real*8 a(*),w(*),u(*),v(*),rv1(*)
        real*8 c,f,g,h,s,x,y,z,scale,anorm
        logical matu,matv
c SUBROUTINE jest przelozona z ALGOLowej wersji SVD z:
c NUM. MATH., 14, 403-420(1970) autorstwa GOLUB i REINSH.
c
c Wyznacza dekompozycje "singular value decomposition" macierzy
c       A = U*S*V    A(m,n)-rzeczywista
c                    U(m,m),S(m,n),V(n,n)
c   UWAGA: iloczyn S(i,i) to wyznacznik macierzy A
c          cond(A)=max{S(i,i)}/min{S(i,i)}=max{w(i)}/min{w(i)}
c Wykorzystuje Hausholder'a bidiagonalizacji algorytm i pewien
c wariant algorytmu QR.
c DANE NA WEJSCIU:
c       a - macierz A(m,n) w postaci wektora zapisana w nim wiersz po wierszu
c               zadeklarowana w programie wywolujacym cz wymiarem najmniej
A(n*m)
c       nm - rozmiar wierszowy macierzy A , a przynajmniej nm>max(m,n)
c       m  - ilosc wierszy macierzy A i U
c       n - ilosc kolumn macierzy A i U i rzad V
c       matu - powinno miec wartosc .TRUE. jesli macierz U w dekompozycji
c               jest oczekiwana, zas .FALSE. w przeciwnym wypadku
c       matv - powinno miec wartosc .TRUE. jesli macierz V w dekompozycji
c               jest oczekiwana, zas .FALSE. w przeciwnym wypadku
c WYNIKI:
c       a - jest niezmieniona ( dopoki nie zapisana przez macierze U lub V)
c       w - zawiera N (nieujemne) osobliwe wartosci A (diagonalne elementy
c               macierzy S). Sa one nieuporzadkowane. Jesli wystapil blad
c               to wartosci osobliwe ierr+1,ierr+2,..,N powinny byc dobre.
c       u - zawiera macierz U (ortogonalnych wektorow - kolumn) z dekompozycji
c               jesli MATU bylo .TRUE. .
c               W przeciwnym wypadku U jest uzywane jako robocza macierz.
c               "U may coincide with A"
c               Jesli blad wystapil to kolumny odpowiadajace prawidlowym
c               wartosciom osobliwym (ierr+1, ierr+2 ...) powinny byc tez
c               prawidlowe.
c               MACIERZ MA ORGANIZACJE WEKTORA ( ZAPISANA W NIM WIERSZ PO
                    WIERSZU)
c  v - zawiera macierz V (ortogonalna) dekompozycji jesli MATV bylo .TRUE.
c               W przeciwnym wypadku jest nie uzywana.
c               "V may also coincide with A"
c               Gdy wykryto blad, kolumny V odpowiadajace poprawnym wartosciom
c               osobliwym powinny byc poprawne.
c               MACIERZ MA ORGANIZACJE WEKTORA ( ZAPISANA W NIM WIERSZ PO
c                 WIERSZU)
c       ierr - wskaznik bledu
c               =0 dla normalnego wykonania procedury
c               =k jesli k-ta osobliwa wartosc nie zostala okreslona po
c               wykonaniu
c                       30 iteracji
c       rv1 - jest macierza robocza
c
c  Pytania i uwagi kierowac do B.S.Garbow, Applied Mathematics Division,
c       Argonne National Labvoratory
c
__________________________________________________________________________
        subroutine uwasi(m,n,a,w,u,v,rv1,uwar,ierr)
c Procedura sprawdza wskaznik uwarunkowania macierzy a zdefiniowany nast:
c       uwar=(max. wart singularna macierzy a)/(min. wart. sing. mac. a)
c macierz a moze byc prostokatna A(m,n) zapisana w postaci wektora wiersz
c       po wierszu.
c Zasada obliczen jest nastepujaca :
c   prowadzi sie rozklad macierzy A = U(m,m)*S(m,n)*V(n,n)
c gdzie U-macierz unormowana o wyznaczniku 1
c       V-macierz unormowana o wyznaczniku 1
c       S-macierz diagonalna o wyznaczniku=det(A)
c         na przekatnej wartosci singularne macierzy A
c UWAGA: jesli m=n , czyli A jest kwadratowa i ktorys z S(i,i)=w(i)=0
c        to macierz A jest osobliwa (det(A)=0)
c       i wtedy ierr=-1
c O ile w trakcie obliczen wykryto blad to ierr > 0
c Macierze U i V nie sa obliczane - traktowane jako robocze
c Macierz W zawiera elementy przekatniowe S, czyli wart singularne A
c Zadane wymiary macierzy a(m*n),w(max(n,m)),rv1(max(n,m)),u(m*m),v(m*n)
c Obliczenia dla Double Precision
        real*8 a(*),w(*),u(*),v(*),rv1(*),uwar
        real*8 wmax,wmin,wa
        integer m,n,nm
        logical matu,matv
__________________________________________________________________________
        subroutine newtset(n,a,f,x1,getfun,getjac,reps,niter,ierr)
c pgm: rozwiazuje nieliniowy układ rownan metoda Newton'a
c       UWAGA: macierze a i f i x1 zmieniaja swa wartosc
c DANE:
c       n - ilosc niewiadomych w układzie rownan ( wraz z mnoznikami
c               Lagrange'a - INTEGER
c       a - macierz Jacobiego dla układu  nieliniowego - tzn macierz pocho-
c               dnych rownan po niewiadomych w punkcie x1 - ma byc wyznacza-
c               na podprogramem getajac(....) - REAL*8
c       f - wektor wartosci residualnych oryginalnego układu  nieliniowego w
c               punkcie x1 - obliczny podprogramem getfun(....) - REAL*8
c       getjac - nazwa podprogramu getjac(....) - EXTERNAL
c       getfun - nazwa podprogramu getfun(....) - EXTERNAL
c       x1 - wektor przyblizenia zerowego rozwiazania na wejsciu - REAL*8
c       reps - mnoznik do maszynowego epsilon dla kryterium zbieznosci
c               kryterium zbieznosci reps*epsmach > norma wekt. przyrostu
c               REAL*8
c       niter - maxymalna ilosc iteracji do wykonania INTEGER
c WYNIKI :
c       x1 - rozwiązanie układu  nieliniowego spelniajacego zadany warunek
c               zbieznosci
c       f - wektor ostatnich residuow
c       reps - norma wektora residuow po zbiegnieciu sie
c       ierr - wskaznik wykonania
c               ierr < 0 byly blady : ierr=-1 macierz Jacobiego osobliwa
c                                     ierr=-2 przekroczenie NITER iteracji
c               ierr > 0 bezblednie : ierr - ilosc wykonanych iteracji
        integer n,ierr,it,niter
        real*8 a(*),f(*),x1(*)
        real*8 deps,eukf,reps,feps
        real*4 eps
__________________________________________________________________________
 

        subroutine epsmach(eps,deps)
c pgm : oblicza tzw epsilon maszynowe tzn najwieksza wartosc dla ktorej
c       jest prawdziwe : 1+eps=1
c       oblicza dla real*4 : eps
c       oblicza dla real*8 : deps
c---------------------------------
        real*8 deps,deps1
        real*4 eps,eps1