text:
!********************************************************************************
!   FORWERR_MULTI: Compare the forward error of FPML, FPML-Quad, and FPML-Comp
!   Author: Thomas R. Cameron, Penn State Erie The Behrend College
!   Last Modified: 31 July 2020
!********************************************************************************
program forwerr_multi
    use fpml, only : main
	use fpml_comp, only: main_comp
	use fpml_quad, only: qp, main_quad
    use mproutines
	implicit none
    ! testing variables
    integer                                     :: deg, j, poly_num
    complex(kind=dp), dimension(:), allocatable :: exact_roots
    ! FPML variables
    integer, parameter                          :: itmax = 30
    integer, allocatable                        :: conv(:)
    real(kind=dp), allocatable                  :: berr(:), cond(:)
	real(kind=qp), allocatable					:: berr_quad(:), cond_quad(:)
    complex(kind=dp), allocatable               :: p(:), roots(:)
	complex(kind=qp), allocatable				:: p_quad(:), roots_quad(:)

    ! Testing: polynomials with multiple and near multiple roots
    open(unit=1,file="data_files/forwerr_multi.dat")
    write(1,'(A)') 'Poly No., FPML, FPML-Quad, FPML-Comp'
    do poly_num=1,7
        write(1, '(I10,A)', advance='no') poly_num, ', '
        ! make polynomial
        call make_poly(poly_num, deg, exact_roots, p)
		! allocate
        allocate(roots(deg), berr(deg), cond(deg), conv(deg))
		allocate(p_quad(deg+1), roots_quad(deg), berr_quad(deg), cond_quad(deg))
        ! FPML
        call main(p, deg, roots, berr, cond, conv, itmax)
        write(1, '(ES15.2, A)', advance='no') max_rel_err(roots,exact_roots,deg), ', '
        ! FPML-Quad
		p_quad = p
        call main_quad(p_quad, deg, roots_quad, berr_quad, cond_quad, conv, itmax)
		roots = roots_quad
        write(1, '(ES15.2, A)', advance='no') max_rel_err(roots,exact_roots,deg), ', '
        ! FPML-Comp
        call main_comp(p, deg, roots, itmax)
        write(1, '(ES15.2, A)') max_rel_err(roots,exact_roots,deg)
        ! deallocate
        deallocate(p, exact_roots, roots, berr, cond, conv)
		deallocate(p_quad, roots_quad, berr_quad, cond_quad)
    end do
    ! close file
    close(1)

contains

		!************************************************
	    !                       make_poly               *
	    !************************************************
	    ! Make special polynomial corresponding to
	    ! poly_num and deg > 1.
	    !************************************************
	    subroutine make_poly(poly_num, deg, exact_roots, p)
	        implicit none
	        ! argument variables
	        integer, intent(in)                         :: poly_num
	        integer, intent(out)                        :: deg
	        complex(kind=dp), allocatable, intent(out)  :: exact_roots(:), p(:)
	        ! local variables
	        integer                                     :: j

	        select case (poly_num)
	            ! case 1
	            case(1)
	            ! Mandelbrot Polynomial, theoretical limiting accuracy 2.53E-5
	            deg = 63
	            allocate(exact_roots(deg), p(deg+1))
	            ! roots
	            exact_roots(1) = cmplx(-0.199909568232701847321062999922E1_dp,0E0_dp,kind=dp)
	            exact_roots(2) = cmplx(-0.199181417254912221573256094986E1_dp,0E0_dp,kind=dp)
	            exact_roots(3) = cmplx(-0.197717958700625738734608852066E1_dp,0E0_dp,kind=dp)
	            exact_roots(4) = cmplx(-0.195370589428439624542762219901E1_dp,0E0_dp,kind=dp)
	            exact_roots(5) = cmplx(-0.192714770936395026246006818895E1_dp,0E0_dp,kind=dp)
	            exact_roots(6) = cmplx(-0.188480357158668179232947809292E1_dp,0E0_dp,kind=dp)
	            exact_roots(7) = cmplx(-0.183231520275122919208489752604E1_dp,0E0_dp,kind=dp)
	            exact_roots(8) = cmplx(-0.176926167027683114607548022864E1_dp,-0.5691950039560031530490018730E-1_dp,kind=dp)
	            exact_roots(9) = cmplx(-0.176926167027683114607548022864E1_dp,0.5691950039560031530490018730E-1_dp,kind=dp)
	            exact_roots(10) = cmplx(-0.167406609147478797156529602917E1_dp,0E0_dp,kind=dp)
	            exact_roots(11) = cmplx(-0.157488913975230096981996555250E1_dp,0E0_dp,kind=dp)
	            exact_roots(12) = cmplx(-0.140844648574007265491757700881E1_dp,-0.13617199730465991568470779361E0_dp,kind=dp)
	            exact_roots(13) = cmplx(-0.140844648574007265491757700881E1_dp,0.13617199730465991568470779361E0_dp,kind=dp)
	            exact_roots(14) = cmplx(-0.1292558061033522087164184706361E1_dp,-0.438198816086631837129737124327E0_dp,kind=dp)
	            exact_roots(15) = cmplx(-0.1292558061033522087164184706361E1_dp,0.438198816086631837129737124327E0_dp,kind=dp)
	            exact_roots(16) = cmplx(-0.1262287281438472543010111941208E1_dp,-0.408104324112690383290160657426E0_dp,kind=dp)
	            exact_roots(17) = cmplx(-0.1262287281438472543010111941208E1_dp,0.408104324112690383290160657426E0_dp,kind=dp)
	            exact_roots(18) = cmplx(-0.1252735884012037946295811002564E1_dp,-0.342470647889750893861875786871E0_dp,kind=dp)
	            exact_roots(19) = cmplx(-0.1252735884012037946295811002564E1_dp,0.342470647889750893861875786871E0_dp,kind=dp)
	            exact_roots(20) = cmplx(-0.1028193852454817599302497455967E1_dp,-0.361376517118561592479460832998E0_dp,kind=dp)
	            exact_roots(21) = cmplx(-0.1028193852454817599302497455967E1_dp,0.361376517118561592479460832998E0_dp,kind=dp)
	            exact_roots(22) = cmplx(-0.623532485956252757990016587001E0_dp,-0.681064414225239608090835812687E0_dp,kind=dp)
	            exact_roots(23) = cmplx(-0.623532485956252757990016587001E0_dp,0.681064414225239608090835812687E0_dp,kind=dp)
	            exact_roots(24) = cmplx(-0.622436295041293587960163506947E0_dp,-0.424878436475629184311574438805E0_dp,kind=dp)
	            exact_roots(25) = cmplx(-0.622436295041293587960163506947E0_dp,0.424878436475629184311574438805E0_dp,kind=dp)
	            exact_roots(26) = cmplx(-0.530827804859942728921477297120E0_dp,-0.668288725559205771444092465565E0_dp,kind=dp)
	            exact_roots(27) = cmplx(-0.530827804859942728921477297120E0_dp,0.668288725559205771444092465565E0_dp,kind=dp)
	            exact_roots(28) = cmplx(-0.272102461488938894219383324518E0_dp,-0.842364690294128145503155708243E0_dp,kind=dp)
	            exact_roots(29) = cmplx(-0.272102461488938894219383324518E0_dp,0.842364690294128145503155708243E0_dp,kind=dp)
	            exact_roots(30) = cmplx(-0.224915951286740054685326255204E0_dp,-0.1116260157454991835001268254245E1_dp,kind=dp)
	            exact_roots(31) = cmplx(-0.224915951286740054685326255204E0_dp,0.1116260157454991835001268254245E1_dp,kind=dp)
	            exact_roots(32) = cmplx(-0.207283835455666412824133850187E0_dp,-0.1117480772494962911373775673122E1_dp,kind=dp)
	            exact_roots(33) = cmplx(-0.207283835455666412824133850187E0_dp,0.1117480772494962911373775673122E1_dp,kind=dp)
	            exact_roots(34) = cmplx(-0.174578221135716969451566432662E0_dp,-0.1071427671454031189229646310220E1_dp,kind=dp)
	            exact_roots(35) = cmplx(-0.174578221135716969451566432662E0_dp,0.1071427671454031189229646310220E1_dp,kind=dp)
	            exact_roots(36) = cmplx(-0.157516053475965356164335109645E0_dp,-0.1109006514113607177971751986155E1_dp,kind=dp)
	            exact_roots(37) = cmplx(-0.157516053475965356164335109645E0_dp,0.1109006514113607177971751986155E1_dp,kind=dp)
	            exact_roots(38) = cmplx(-0.127499973546363000199539565346E0_dp,-0.987460909489456792207407680793E0_dp,kind=dp)
	            exact_roots(39) = cmplx(-0.127499973546363000199539565346E0_dp,0.987460909489456792207407680793E0_dp,kind=dp)
	            exact_roots(40) = cmplx(-0.14233481920354066767761845314E-1_dp,-0.1032914775213644109395013402655E1_dp,kind=dp)
	            exact_roots(41) = cmplx(-0.14233481920354066767761845314E-1_dp,0.1032914775213644109395013402655E1_dp,kind=dp)
	            exact_roots(42) = cmplx(-0.6983568496261391817966491075E-2_dp,-0.1003603862288289548530704966951E1_dp,kind=dp)
	            exact_roots(43) = cmplx(-0.6983568496261391817966491075E-2_dp,0.1003603862288289548530704966951E1_dp,kind=dp)
	            exact_roots(44) = cmplx(0.14895466603687646529815779209E-1_dp,-0.848148761908416527719331111783E0_dp,kind=dp)
	            exact_roots(45) = cmplx(0.14895466603687646529815779209E-1_dp,0.848148761908416527719331111783E0_dp,kind=dp)
	            exact_roots(46) = cmplx(0.121192786105906486314704443411E0_dp,-0.610611692210754211675387244150E0_dp,kind=dp)
	            exact_roots(47) = cmplx(0.121192786105906486314704443411E0_dp,0.610611692210754211675387244150E0_dp,kind=dp)
	            exact_roots(48) = cmplx(0.352482539722363278193253964052E0_dp,-0.698337239583330331258141954760E0_dp,kind=dp)
	            exact_roots(49) = cmplx(0.352482539722363278193253964052E0_dp,0.698337239583330331258141954760E0_dp,kind=dp)
	            exact_roots(50) = cmplx(0.376008681846767559704804317729E0_dp,-0.144749371321632864747110182013E0_dp,kind=dp)
	            exact_roots(51) = cmplx(0.376008681846767559704804317729E0_dp,0.144749371321632864747110182013E0_dp,kind=dp)
	            exact_roots(52) = cmplx(0.376893240379311323690004017969E0_dp,-0.678568693190448141957540792997E0_dp,kind=dp)
	            exact_roots(53) = cmplx(0.376893240379311323690004017969E0_dp,0.678568693190448141957540792997E0_dp,kind=dp)
	            exact_roots(54) = cmplx(0.386539176596158026508293086904E0_dp,-0.569324711303102903213792357135E0_dp,kind=dp)
	            exact_roots(55) = cmplx(0.386539176596158026508293086904E0_dp,0.569324711303102903213792357135E0_dp,kind=dp)
	            exact_roots(56) = cmplx(0.412916024722700479197334566383E0_dp,-0.614806760143385694954549720401E0_dp,kind=dp)
	            exact_roots(57) = cmplx(0.412916024722700479197334566383E0_dp,0.614806760143385694954549720401E0_dp,kind=dp)
	            exact_roots(58) = cmplx(0.432376192641994507824669648087E0_dp,-0.226759904435348618697876559972E0_dp,kind=dp)
	            exact_roots(59) = cmplx(0.432376192641994507824669648087E0_dp,0.226759904435348618697876559972E0_dp,kind=dp)
	            exact_roots(60) = cmplx(0.452774498724915493508803077733E0_dp,-0.396170128033165002412596877271E0_dp,kind=dp)
	            exact_roots(61) = cmplx(0.452774498724915493508803077733E0_dp,0.396170128033165002412596877271E0_dp,kind=dp)
	            exact_roots(62) = cmplx(0.456823285823316651283953236253E0_dp,-0.347758700883481983632188723200E0_dp,kind=dp)
	            exact_roots(63) = cmplx(0.456823285823316651283953236253E0_dp,0.347758700883481983632188723200E0_dp,kind=dp)
	            ! polynomial
	            p(1) = cmplx(1E0_dp,0E0_dp,kind=dp)
	            p(2) = cmplx(1E0_dp,0E0_dp,kind=dp)
	            p(3) = cmplx(2E0_dp,0E0_dp,kind=dp)
	            p(4) = cmplx(5E0_dp,0E0_dp,kind=dp)
	            p(5) = cmplx(14E0_dp,0E0_dp,kind=dp)
	            p(6) = cmplx(42E0_dp,0E0_dp,kind=dp)
	            p(7) = cmplx(132E0_dp,0E0_dp,kind=dp)
	            p(8) = cmplx(365E0_dp,0E0_dp,kind=dp)
	            p(9) = cmplx(950E0_dp,0E0_dp,kind=dp)
	            p(10) = cmplx(2398E0_dp,0E0_dp,kind=dp)
	            p(11) = cmplx(5916E0_dp,0E0_dp,kind=dp)
	            p(12) = cmplx(14290E0_dp,0E0_dp,kind=dp)
	            p(13) = cmplx(33708E0_dp,0E0_dp,kind=dp)
	            p(14) = cmplx(77684E0_dp,0E0_dp,kind=dp)
	            p(15) = cmplx(175048E0_dp,0E0_dp,kind=dp)
	            p(16) = cmplx(385741E0_dp,0E0_dp,kind=dp)
	            p(17) = cmplx(831014E0_dp,0E0_dp,kind=dp)
	            p(18) = cmplx(1749654E0_dp,0E0_dp,kind=dp)
	            p(19) = cmplx(3598964E0_dp,0E0_dp,kind=dp)
	            p(20) = cmplx(7228014E0_dp,0E0_dp,kind=dp)
	            p(21) = cmplx(14162220E0_dp,0E0_dp,kind=dp)
	            p(22) = cmplx(27049196E0_dp,0E0_dp,kind=dp)
	            p(23) = cmplx(50323496E0_dp,0E0_dp,kind=dp)
	            p(24) = cmplx(91143114E0_dp,0E0_dp,kind=dp)
	            p(25) = cmplx(160617860E0_dp,0E0_dp,kind=dp)
	            p(26) = cmplx(275276716E0_dp,0E0_dp,kind=dp)
	            p(27) = cmplx(458591432E0_dp,0E0_dp,kind=dp)
	            p(28) = cmplx(742179284E0_dp,0E0_dp,kind=dp)
	            p(29) = cmplx(1166067016E0_dp,0E0_dp,kind=dp)
	            p(30) = cmplx(1777171560E0_dp,0E0_dp,kind=dp)
	            p(31) = cmplx(2625062128E0_dp,0E0_dp,kind=dp)
	            p(32) = cmplx(3754272037E0_dp,0E0_dp,kind=dp)
	            p(33) = cmplx(5193067630E0_dp,0E0_dp,kind=dp)
	            p(34) = cmplx(6939692682E0_dp,0E0_dp,kind=dp)
	            p(35) = cmplx(8948546308E0_dp,0E0_dp,kind=dp)
	            p(36) = cmplx(11120136162E0_dp,0E0_dp,kind=dp)
	            p(37) = cmplx(13299362332E0_dp,0E0_dp,kind=dp)
	            p(38) = cmplx(15286065700E0_dp,0E0_dp,kind=dp)
	            p(39) = cmplx(16859410792E0_dp,0E0_dp,kind=dp)
	            p(40) = cmplx(17813777994E0_dp,0E0_dp,kind=dp)
	            p(41) = cmplx(17999433372E0_dp,0E0_dp,kind=dp)
	            p(42) = cmplx(17357937708E0_dp,0E0_dp,kind=dp)
	            p(43) = cmplx(15941684776E0_dp,0E0_dp,kind=dp)
	            p(44) = cmplx(13910043524E0_dp,0E0_dp,kind=dp)
	            p(45) = cmplx(11500901864E0_dp,0E0_dp,kind=dp)
	            p(46) = cmplx(8984070856E0_dp,0E0_dp,kind=dp)
	            p(47) = cmplx(6609143792E0_dp,0E0_dp,kind=dp)
	            p(48) = cmplx(4562339774E0_dp,0E0_dp,kind=dp)
	            p(49) = cmplx(2943492972E0_dp,0E0_dp,kind=dp)
	            p(50) = cmplx(1766948340E0_dp,0E0_dp,kind=dp)
	            p(51) = cmplx(981900168E0_dp,0E0_dp,kind=dp)
	            p(52) = cmplx(502196500E0_dp,0E0_dp,kind=dp)
	            p(53) = cmplx(234813592E0_dp,0E0_dp,kind=dp)
	            p(54) = cmplx(99582920E0_dp,0E0_dp,kind=dp)
	            p(55) = cmplx(37945904E0_dp,0E0_dp,kind=dp)
	            p(56) = cmplx(12843980E0_dp,0E0_dp,kind=dp)
	            p(57) = cmplx(3807704E0_dp,0E0_dp,kind=dp)
	            p(58) = cmplx(971272E0_dp,0E0_dp,kind=dp)
	            p(59) = cmplx(208336E0_dp,0E0_dp,kind=dp)
	            p(60) = cmplx(36440E0_dp,0E0_dp,kind=dp)
	            p(61) = cmplx(4976E0_dp,0E0_dp,kind=dp)
	            p(62) = cmplx(496E0_dp,0E0_dp,kind=dp)
	            p(63) = cmplx(32E0_dp,0E0_dp,kind=dp)
	            p(64) = cmplx(1E0_dp,0E0_dp,kind=dp)
	            ! case 2
	            case(2)
	            ! Kameny Polynomial c=10, theoretical limiting accuracy 1.11E-16
	            deg = 9
	            allocate(exact_roots(deg), p(deg+1))
	            ! roots
	            exact_roots(1) = cmplx(-0.250708394126194522431084083966E1_dp,0E0_dp,kind=dp)
	            exact_roots(2) = cmplx(-0.77475933255040504312683246509E0_dp,-0.239347474286129783412456829950E1_dp,kind=dp)
	            exact_roots(3) = cmplx(-0.77475933255040504312683246509E0_dp,0.239347474286129783412456829950E1_dp,kind=dp)
	            exact_roots(4) = cmplx(-0.173313478195944654943157024747E0_dp,0E0_dp,kind=dp)
	            exact_roots(5) = cmplx(-0.173097223325097844852013748504E0_dp,0E0_dp,kind=dp)
	            exact_roots(6) = cmplx(0.173204810760521109164135307321E0_dp,-0.108125591879847889809932772E-3_dp,kind=dp)
	            exact_roots(7) = cmplx(0.173204810760521109164135307321E0_dp,0.108125591879847889809932772E-3_dp,kind=dp)
	            exact_roots(8) = cmplx(0.202830184318137779601570296422E1_dp,-0.147928157900441456648083759833E1_dp,kind=dp)
	            exact_roots(9) = cmplx(0.202830184318137779601570296422E1_dp,0.147928157900441456648083759833E1_dp,kind=dp)
	            ! polynomial
	            p(1) = cmplx(9E0_dp,0E0_dp,kind=dp)
	            p(2) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(3) = cmplx(-600E0_dp,0E0_dp,kind=dp)
	            p(4) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(5) = cmplx(10000E0_dp,0E0_dp,kind=dp)
	            p(6) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(7) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(8) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(9) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(10) = cmplx(100E0_dp,0E0_dp,kind=dp)
	            ! case 3
	            case(3)
	            ! Kameny Polynomial c=10^3, theoretical limiting accuracy 1.11E-16
	            deg = 9
	            allocate(exact_roots(deg), p(deg+1))
	            ! roots
	            exact_roots(1) = cmplx(-0.158489318488962525196907753240E2_dp,0E0_dp,kind=dp)
	            exact_roots(2) = cmplx(-0.48975892839992994551324676957E1_dp,-0.150732300552288401309352587876E2_dp,kind=dp)
	            exact_roots(3) = cmplx(-0.48975892839992994551324676957E1_dp,0.150732300552288401309352587876E2_dp,kind=dp)
	            exact_roots(4) = cmplx(-0.173205080767700380719051028416E-2_dp,0E0_dp,kind=dp)
	            exact_roots(5) = cmplx(-0.173205080746075077991838239885E-2_dp,0E0_dp,kind=dp)
	            exact_roots(6) = cmplx(0.173205080756887729350044634151E-2_dp,-0.10812651363606394264E-12_dp,kind=dp)
	            exact_roots(7) = cmplx(0.173205080756887729350044634151E-2_dp,0.10812651363606394264E-12_dp,kind=dp)
	            exact_roots(8) = cmplx(0.128220552084474257149779093577E2_dp,-0.93157684943778791903260859371E1_dp,kind=dp)
	            exact_roots(9) = cmplx(0.128220552084474257149779093577E2_dp,0.93157684943778791903260859371E1_dp,kind=dp)
	            ! polynomial
	            p(1) = cmplx(9E0_dp,0E0_dp,kind=dp)
	            p(2) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(3) = cmplx(-6000000E0_dp,0E0_dp,kind=dp)
	            p(4) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(5) = cmplx(1000000000000E0_dp,0E0_dp,kind=dp)
	            p(6) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(7) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(8) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(9) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(10) = cmplx(1000000E0_dp,0E0_dp,kind=dp)
	            ! case 4
	            case(4)
	            ! Kameny Polynomial c=10^6, theoretical limiting accuracy 8.09E-9
	            deg = 9
	            allocate(exact_roots(deg), p(deg+1))
	            ! roots
	            exact_roots(1) = cmplx(-0.251188643150958006331217160138E3_dp,0E0_dp,kind=dp)
	            exact_roots(2) = cmplx(-0.77621559527630266049907587459E2_dp,-0.238894595888056617951698213041E3_dp,kind=dp)
	            exact_roots(3) = cmplx(-0.77621559527630266049907587459E2_dp,0.238894595888056617951698213041E3_dp,kind=dp)
	            exact_roots(4) = cmplx(-0.173205080756887729353086560209E-5_dp,0E0_dp,kind=dp)
	            exact_roots(5) = cmplx(-0.173205080756887729352402708092E-5_dp,0E0_dp,kind=dp)
	            exact_roots(6) = cmplx(0.173205080756887729352744634151E-5_dp,-0.341926059E-26_dp,kind=dp)
	            exact_roots(7) = cmplx(0.173205080756887729352744634151E-5_dp,0.341926059E-26_dp,kind=dp)
	            exact_roots(8) = cmplx(0.203215881103109269215516167528E3_dp,-0.147644979987489859882850769633E3_dp,kind=dp)
	            exact_roots(9) = cmplx(0.203215881103109269215516167528E3_dp,0.147644979987489859882850769633E3_dp,kind=dp)
	            ! polynomial
	            p(1) = cmplx(9E0_dp,0E0_dp,kind=dp)
	            p(2) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(3) = cmplx(-6000000000000E0_dp,0E0_dp,kind=dp)
	            p(4) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(5) = cmplx(1000000000000000000000000E0_dp,0E0_dp,kind=dp)
	            p(6) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(7) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(8) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(9) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(10) = cmplx(1000000000000E0_dp,0E0_dp,kind=dp)
	            ! case 5
	            case(5)
	            ! Multiplicity I
	            deg = 55
	            allocate(exact_roots(deg), p(deg+1))
	            ! roots
	            exact_roots(1) = cmplx(-0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(2) = cmplx(-0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(3) = cmplx(-0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(4) = cmplx(-0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(5) = cmplx(-0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(6) = cmplx(-0.951677511788434760780566076462E0_dp,-0.156577036126765257092467586019E0_dp,kind=dp)
	            exact_roots(7) = cmplx(-0.951677511788434760780566076462E0_dp,0.156577036126765257092467586019E0_dp,kind=dp)
	            exact_roots(8) = cmplx(-0.947111794804360643443628545664E0_dp,-0.45987543603302557306898551016E-1_dp,kind=dp)
	            exact_roots(9) = cmplx(-0.947111794804360643443628545664E0_dp,0.45987543603302557306898551016E-1_dp,kind=dp)
	            exact_roots(10) = cmplx(-0.935147084823051581325738268702E0_dp,-0.276366395363089095701385358741E0_dp,kind=dp)
	            exact_roots(11) = cmplx(-0.935147084823051581325738268702E0_dp,0.276366395363089095701385358741E0_dp,kind=dp)
	            exact_roots(12) = cmplx(-0.899451950794801943937381293773E0_dp,-0.394568367012592113985865567971E0_dp,kind=dp)
	            exact_roots(13) = cmplx(-0.899451950794801943937381293773E0_dp,0.394568367012592113985865567971E0_dp,kind=dp)
	            exact_roots(14) = cmplx(-0.846928410223891427281972003960E0_dp,-0.507591758019284469647521636947E0_dp,kind=dp)
	            exact_roots(15) = cmplx(-0.846928410223891427281972003960E0_dp,0.507591758019284469647521636947E0_dp,kind=dp)
	            exact_roots(16) = cmplx(-0.779195724062158658036386916327E0_dp,-0.613081457400403535591975580892E0_dp,kind=dp)
	            exact_roots(17) = cmplx(-0.779195724062158658036386916327E0_dp,0.613081457400403535591975580892E0_dp,kind=dp)
	            exact_roots(18) = cmplx(-0.697743392981121831801589116937E0_dp,-0.709081007070975090262670227767E0_dp,kind=dp)
	            exact_roots(19) = cmplx(-0.697743392981121831801589116937E0_dp,0.709081007070975090262670227767E0_dp,kind=dp)
	            exact_roots(20) = cmplx(-0.604121886825094827201080735988E0_dp,-0.793887314714735149770017023538E0_dp,kind=dp)
	            exact_roots(21) = cmplx(-0.604121886825094827201080735988E0_dp,0.793887314714735149770017023538E0_dp,kind=dp)
	            exact_roots(22) = cmplx(-0.500000000000000000000000000000E0_dp,-0.866025403784438646763723170753E0_dp,kind=dp)
	            exact_roots(23) = cmplx(-0.500000000000000000000000000000E0_dp,0.866025403784438646763723170753E0_dp,kind=dp)
	            exact_roots(24) = cmplx(-0.387174132037879582862879778205E0_dp,-0.924251993602053902039550620488E0_dp,kind=dp)
	            exact_roots(25) = cmplx(-0.387174132037879582862879778205E0_dp,0.924251993602053902039550620488E0_dp,kind=dp)
	            exact_roots(26) = cmplx(-0.267556959428167474972897159502E0_dp,-0.967566373494130224905327435791E0_dp,kind=dp)
	            exact_roots(27) = cmplx(-0.267556959428167474972897159502E0_dp,0.967566373494130224905327435791E0_dp,kind=dp)
	            exact_roots(28) = cmplx(-0.143155839532593341810009310219E0_dp,-0.995221832724034632237407792227E0_dp,kind=dp)
	            exact_roots(29) = cmplx(-0.143155839532593341810009310219E0_dp,0.995221832724034632237407792227E0_dp,kind=dp)
	            exact_roots(30) = cmplx(-0.16045517947534302681825863292E-1_dp,-0.1006735041626386010701028433184E1_dp,kind=dp)
	            exact_roots(31) = cmplx(-0.16045517947534302681825863292E-1_dp,0.1006735041626386010701028433184E1_dp,kind=dp)
	            exact_roots(32) = cmplx(0.111662481142841735790134563805E0_dp,-0.1001892206008032435299692737060E1_dp,kind=dp)
	            exact_roots(33) = cmplx(0.111662481142841735790134563805E0_dp,0.1001892206008032435299692737060E1_dp,kind=dp)
	            exact_roots(34) = cmplx(0.237852321567861477202899799722E0_dp,-0.980751410679651646035316169672E0_dp,kind=dp)
	            exact_roots(35) = cmplx(0.237852321567861477202899799722E0_dp,0.980751410679651646035316169672E0_dp,kind=dp)
	            exact_roots(36) = cmplx(0.360437139274099643543556447604E0_dp,-0.943640865401314173256925688139E0_dp,kind=dp)
	            exact_roots(37) = cmplx(0.360437139274099643543556447604E0_dp,0.943640865401314173256925688139E0_dp,kind=dp)
	            exact_roots(38) = cmplx(0.477392314844102847110474634171E0_dp,-0.891152945345826129924640822366E0_dp,kind=dp)
	            exact_roots(39) = cmplx(0.477392314844102847110474634171E0_dp,0.891152945345826129924640822366E0_dp,kind=dp)
	            exact_roots(40) = cmplx(0.586788023132224235580014756671E0_dp,-0.824134044230047196220784957068E0_dp,kind=dp)
	            exact_roots(41) = cmplx(0.586788023132224235580014756671E0_dp,0.824134044230047196220784957068E0_dp,kind=dp)
	            exact_roots(42) = cmplx(0.686820425487486753027867016412E0_dp,-0.743670357958749198637378689329E0_dp,kind=dp)
	            exact_roots(43) = cmplx(0.686820425487486753027867016412E0_dp,0.743670357958749198637378689329E0_dp,kind=dp)
	            exact_roots(44) = cmplx(0.775840929682842739762490529077E0_dp,-0.651069801113127816898897251309E0_dp,kind=dp)
	            exact_roots(45) = cmplx(0.775840929682842739762490529077E0_dp,0.651069801113127816898897251309E0_dp,kind=dp)
	            exact_roots(46) = cmplx(0.852382996720766159190958138062E0_dp,-0.547840332639938504859225568105E0_dp,kind=dp)
	            exact_roots(47) = cmplx(0.852382996720766159190958138062E0_dp,0.547840332639938504859225568105E0_dp,kind=dp)
	            exact_roots(48) = cmplx(0.915186027756796108167454007547E0_dp,-0.435665032553729539118325217757E0_dp,kind=dp)
	            exact_roots(49) = cmplx(0.915186027756796108167454007547E0_dp,0.435665032553729539118325217757E0_dp,kind=dp)
	            exact_roots(50) = cmplx(0.963215921119370289466626861950E0_dp,-0.316374328936911136376493374596E0_dp,kind=dp)
	            exact_roots(51) = cmplx(0.963215921119370289466626861950E0_dp,0.316374328936911136376493374596E0_dp,kind=dp)
	            exact_roots(52) = cmplx(0.995681949864471509895340077408E0_dp,-0.191915823987060561907739223463E0_dp,kind=dp)
	            exact_roots(53) = cmplx(0.995681949864471509895340077408E0_dp,0.191915823987060561907739223463E0_dp,kind=dp)
	            exact_roots(54) = cmplx(0.1012049674656226877398138236602E1_dp,-0.64322209034046437066803422584E-1_dp,kind=dp)
	            exact_roots(55) = cmplx(0.1012049674656226877398138236602E1_dp,0.64322209034046437066803422584E-1_dp,kind=dp)
	            ! polynomial
	            p(1) = cmplx(1E0_dp,0E0_dp,kind=dp)
	            p(2) = cmplx(6E0_dp,0E0_dp,kind=dp)
	            p(3) = cmplx(15E0_dp,0E0_dp,kind=dp)
	            p(4) = cmplx(20E0_dp,0E0_dp,kind=dp)
	            p(5) = cmplx(15E0_dp,0E0_dp,kind=dp)
	            p(6) = cmplx(6E0_dp,0E0_dp,kind=dp)
	            p(7) = cmplx(1E0_dp,0E0_dp,kind=dp)
	            p(8) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(9) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(10) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(11) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(12) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(13) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(14) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(15) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(16) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(17) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(18) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(19) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(20) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(21) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(22) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(23) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(24) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(25) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(26) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(27) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(28) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(29) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(30) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(31) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(32) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(33) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(34) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(35) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(36) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(37) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(38) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(39) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(40) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(41) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(42) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(43) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(44) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(45) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(46) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(47) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(48) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(49) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(50) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(51) = cmplx(1E0_dp,0E0_dp,kind=dp)
	            p(52) = cmplx(5E0_dp,0E0_dp,kind=dp)
	            p(53) = cmplx(10E0_dp,0E0_dp,kind=dp)
	            p(54) = cmplx(10E0_dp,0E0_dp,kind=dp)
	            p(55) = cmplx(5E0_dp,0E0_dp,kind=dp)
	            p(56) = cmplx(1E0_dp,0E0_dp,kind=dp)
	            ! case 6
	            case(6)
	            ! Multiplicity II
	            deg = 62
	            allocate(exact_roots(deg), p(deg+1))
	            ! roots
	            exact_roots(1) = cmplx(-0.998026728428271561952336806863E0_dp,-0.62790519529313376076178224566E-1_dp,kind=dp)
	            exact_roots(2) = cmplx(-0.998026728428271561952336806863E0_dp,0.62790519529313376076178224566E-1_dp,kind=dp)
	            exact_roots(3) = cmplx(-0.982287250728688681085641742865E0_dp,-0.187381314585724630542550734447E0_dp,kind=dp)
	            exact_roots(4) = cmplx(-0.982287250728688681085641742865E0_dp,0.187381314585724630542550734447E0_dp,kind=dp)
	            exact_roots(5) = cmplx(-0.951056516295153572116439333379E0_dp,-0.309016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(6) = cmplx(-0.951056516295153572116439333379E0_dp,0.309016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(7) = cmplx(-0.904827052466019527713668647933E0_dp,-0.425779291565072648862502445744E0_dp,kind=dp)
	            exact_roots(8) = cmplx(-0.904827052466019527713668647933E0_dp,0.425779291565072648862502445744E0_dp,kind=dp)
	            exact_roots(9) = cmplx(-0.844327925502015078548558063967E0_dp,-0.535826794978996618271308767868E0_dp,kind=dp)
	            exact_roots(10) = cmplx(-0.844327925502015078548558063967E0_dp,0.535826794978996618271308767868E0_dp,kind=dp)
	            exact_roots(11) = cmplx(-0.770513242775789230803009636396E0_dp,-0.637423989748689710176712811676E0_dp,kind=dp)
	            exact_roots(12) = cmplx(-0.770513242775789230803009636396E0_dp,0.637423989748689710176712811676E0_dp,kind=dp)
	            exact_roots(13) = cmplx(-0.684547105928688673732283357621E0_dp,-0.728968627421411523146730319055E0_dp,kind=dp)
	            exact_roots(14) = cmplx(-0.684547105928688673732283357621E0_dp,0.728968627421411523146730319055E0_dp,kind=dp)
	            exact_roots(15) = cmplx(-0.587785252292473129168705954639E0_dp,-0.809016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(16) = cmplx(-0.587785252292473129168705954639E0_dp,0.809016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(17) = cmplx(-0.50000000000000000000000000000E0_dp,-0.217944947177033677611849099193E1_dp,kind=dp)
	            exact_roots(18) = cmplx(-0.50000000000000000000000000000E0_dp,-0.217944947177033677611849099193E1_dp,kind=dp)
	            exact_roots(19) = cmplx(-0.50000000000000000000000000000E0_dp,-0.217944947177033677611849099193E1_dp,kind=dp)
	            exact_roots(20) = cmplx(-0.50000000000000000000000000000E0_dp,0.217944947177033677611849099193E1_dp,kind=dp)
	            exact_roots(21) = cmplx(-0.50000000000000000000000000000E0_dp,0.217944947177033677611849099193E1_dp,kind=dp)
	            exact_roots(22) = cmplx(-0.50000000000000000000000000000E0_dp,0.217944947177033677611849099193E1_dp,kind=dp)
	            exact_roots(23) = cmplx(-0.481753674101715274987191502872E0_dp,-0.876306680043863587308115903922E0_dp,kind=dp)
	            exact_roots(24) = cmplx(-0.481753674101715274987191502872E0_dp,0.876306680043863587308115903922E0_dp,kind=dp)
	            exact_roots(25) = cmplx(-0.368124552684677959156947147493E0_dp,-0.929776485888251403660942556222E0_dp,kind=dp)
	            exact_roots(26) = cmplx(-0.368124552684677959156947147493E0_dp,0.929776485888251403660942556222E0_dp,kind=dp)
	            exact_roots(27) = cmplx(-0.248689887164854788242283746006E0_dp,-0.968583161128631119490168375465E0_dp,kind=dp)
	            exact_roots(28) = cmplx(-0.248689887164854788242283746006E0_dp,0.968583161128631119490168375465E0_dp,kind=dp)
	            exact_roots(29) = cmplx(-0.125333233564304245373118759817E0_dp,-0.992114701314477831049793042786E0_dp,kind=dp)
	            exact_roots(30) = cmplx(-0.125333233564304245373118759817E0_dp,0.992114701314477831049793042786E0_dp,kind=dp)
	            exact_roots(31) = cmplx(0E0_dp,-0.100000000000000000000000000000E1_dp,kind=dp)
	            exact_roots(32) = cmplx(0E0_dp,0.100000000000000000000000000000E1_dp,kind=dp)
	            exact_roots(33) = cmplx(0.125333233564304245373118759817E0_dp,-0.992114701314477831049793042786E0_dp,kind=dp)
	            exact_roots(34) = cmplx(0.125333233564304245373118759817E0_dp,0.992114701314477831049793042786E0_dp,kind=dp)
	            exact_roots(35) = cmplx(0.248689887164854788242283746006E0_dp,-0.968583161128631119490168375465E0_dp,kind=dp)
	            exact_roots(36) = cmplx(0.248689887164854788242283746006E0_dp,0.968583161128631119490168375465E0_dp,kind=dp)
	            exact_roots(37) = cmplx(0.333333333333333333333333333333E0_dp,0E0_dp,kind=dp)
	            exact_roots(38) = cmplx(0.333333333333333333333333333333E0_dp,0E0_dp,kind=dp)
	            exact_roots(39) = cmplx(0.368124552684677959156947147493E0_dp,-0.929776485888251403660942556222E0_dp,kind=dp)
	            exact_roots(40) = cmplx(0.368124552684677959156947147493E0_dp,0.929776485888251403660942556222E0_dp,kind=dp)
	            exact_roots(41) = cmplx(0.481753674101715274987191502872E0_dp,-0.876306680043863587308115903922E0_dp,kind=dp)
	            exact_roots(42) = cmplx(0.481753674101715274987191502872E0_dp,0.876306680043863587308115903922E0_dp,kind=dp)
	            exact_roots(43) = cmplx(0.587785252292473129168705954639E0_dp,-0.809016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(44) = cmplx(0.587785252292473129168705954639E0_dp,0.809016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(45) = cmplx(0.684547105928688673732283357621E0_dp,-0.728968627421411523146730319055E0_dp,kind=dp)
	            exact_roots(46) = cmplx(0.684547105928688673732283357621E0_dp,0.728968627421411523146730319055E0_dp,kind=dp)
	            exact_roots(47) = cmplx(0.770513242775789230803009636396E0_dp,-0.637423989748689710176712811676E0_dp,kind=dp)
	            exact_roots(48) = cmplx(0.770513242775789230803009636396E0_dp,0.637423989748689710176712811676E0_dp,kind=dp)
	            exact_roots(49) = cmplx(0.844327925502015078548558063967E0_dp,-0.535826794978996618271308767868E0_dp,kind=dp)
	            exact_roots(50) = cmplx(0.844327925502015078548558063967E0_dp,0.535826794978996618271308767868E0_dp,kind=dp)
	            exact_roots(51) = cmplx(0.904827052466019527713668647933E0_dp,-0.425779291565072648862502445744E0_dp,kind=dp)
	            exact_roots(52) = cmplx(0.904827052466019527713668647933E0_dp,0.425779291565072648862502445744E0_dp,kind=dp)
	            exact_roots(53) = cmplx(0.951056516295153572116439333379E0_dp,-0.309016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(54) = cmplx(0.951056516295153572116439333379E0_dp,0.309016994374947424102293417183E0_dp,kind=dp)
	            exact_roots(55) = cmplx(0.982287250728688681085641742865E0_dp,-0.187381314585724630542550734447E0_dp,kind=dp)
	            exact_roots(56) = cmplx(0.982287250728688681085641742865E0_dp,0.187381314585724630542550734447E0_dp,kind=dp)
	            exact_roots(57) = cmplx(0.998026728428271561952336806863E0_dp,-0.62790519529313376076178224566E-1_dp,kind=dp)
	            exact_roots(58) = cmplx(0.998026728428271561952336806863E0_dp,0.62790519529313376076178224566E-1_dp,kind=dp)
	            exact_roots(59) = cmplx(0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(60) = cmplx(0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(61) = cmplx(0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(62) = cmplx(0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            ! polynomial
	            p(1) = cmplx(125E0_dp,0E0_dp,kind=dp)
	            p(2) = cmplx(-1175E0_dp,0E0_dp,kind=dp)
	            p(3) = cmplx(4215E0_dp,0E0_dp,kind=dp)
	            p(4) = cmplx(-7444E0_dp,0E0_dp,kind=dp)
	            p(5) = cmplx(7393E0_dp,0E0_dp,kind=dp)
	            p(6) = cmplx(-5133E0_dp,0E0_dp,kind=dp)
	            p(7) = cmplx(3402E0_dp,0E0_dp,kind=dp)
	            p(8) = cmplx(-1917E0_dp,0E0_dp,kind=dp)
	            p(9) = cmplx(741E0_dp,0E0_dp,kind=dp)
	            p(10) = cmplx(-316E0_dp,0E0_dp,kind=dp)
	            p(11) = cmplx(115E0_dp,0E0_dp,kind=dp)
	            p(12) = cmplx(-15E0_dp,0E0_dp,kind=dp)
	            p(13) = cmplx(9E0_dp,0E0_dp,kind=dp)
	            p(14) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(15) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(16) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(17) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(18) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(19) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(20) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(21) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(22) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(23) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(24) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(25) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(26) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(27) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(28) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(29) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(30) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(31) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(32) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(33) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(34) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(35) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(36) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(37) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(38) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(39) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(40) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(41) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(42) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(43) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(44) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(45) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(46) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(47) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(48) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(49) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(50) = cmplx(0E0_dp,0E0_dp,kind=dp)
	            p(51) = cmplx(125E0_dp,0E0_dp,kind=dp)
	            p(52) = cmplx(-1175E0_dp,0E0_dp,kind=dp)
	            p(53) = cmplx(4215E0_dp,0E0_dp,kind=dp)
	            p(54) = cmplx(-7444E0_dp,0E0_dp,kind=dp)
	            p(55) = cmplx(7393E0_dp,0E0_dp,kind=dp)
	            p(56) = cmplx(-5133E0_dp,0E0_dp,kind=dp)
	            p(57) = cmplx(3402E0_dp,0E0_dp,kind=dp)
	            p(58) = cmplx(-1917E0_dp,0E0_dp,kind=dp)
	            p(59) = cmplx(741E0_dp,0E0_dp,kind=dp)
	            p(60) = cmplx(-316E0_dp,0E0_dp,kind=dp)
	            p(61) = cmplx(115E0_dp,0E0_dp,kind=dp)
	            p(62) = cmplx(-15E0_dp,0E0_dp,kind=dp)
	            p(63) = cmplx(9E0_dp,0E0_dp,kind=dp)
	            ! case 7
	            case(7)
	            ! Multiplicity III
	            deg = 17
	            allocate(exact_roots(deg), p(deg+1))
	            ! roots
	            exact_roots(1) = cmplx(0.100000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(2) = cmplx(0.200000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(3) = cmplx(0.300000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(4) = cmplx(0.400000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(5) = cmplx(0.500000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(6) = cmplx(0.600000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(7) = cmplx(0.700000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(8) = cmplx(0.800000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(9) = cmplx(0.900000000000000000000000000000E1_dp,0E0_dp,kind=dp)
	            exact_roots(10) = cmplx(0.100000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            exact_roots(11) = cmplx(0.110000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            exact_roots(12) = cmplx(0.120000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            exact_roots(13) = cmplx(0.130000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            exact_roots(14) = cmplx(0.140000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            exact_roots(15) = cmplx(0.150000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            exact_roots(16) = cmplx(0.150000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            exact_roots(17) = cmplx(0.150000000000000000000000000000E2_dp,0E0_dp,kind=dp)
	            ! polynomial
	            p(1) = cmplx(-294226732800000E0_dp,0E0_dp,kind=dp)
	            p(2) = cmplx(1015541906400000E0_dp,0E0_dp,kind=dp)
	            p(3) = cmplx(-1518791527728000E0_dp,0E0_dp,kind=dp)
	            p(4) = cmplx(1327137724803600E0_dp,0E0_dp,kind=dp)
	            p(5) = cmplx(-766908691489440E0_dp,0E0_dp,kind=dp)
	            p(6) = cmplx(313437620164824E0_dp,0E0_dp,kind=dp)
	            p(7) = cmplx(-94377698961000E0_dp,0E0_dp,kind=dp)
	            p(8) = cmplx(21485772576905E0_dp,0E0_dp,kind=dp)
	            p(9) = cmplx(-3758453397270E0_dp,0E0_dp,kind=dp)
	            p(10) = cmplx(509681511053E0_dp,0E0_dp,kind=dp)
	            p(11) = cmplx(-53726158200E0_dp,0E0_dp,kind=dp)
	            p(12) = cmplx(4387265090E0_dp,0E0_dp,kind=dp)
	            p(13) = cmplx(-274687140E0_dp,0E0_dp,kind=dp)
	            p(14) = cmplx(12932122E0_dp,0E0_dp,kind=dp)
	            p(15) = cmplx(-442800E0_dp,0E0_dp,kind=dp)
	            p(16) = cmplx(10405E0_dp,0E0_dp,kind=dp)
	            p(17) = cmplx(-150E0_dp,0E0_dp,kind=dp)
	            p(18) = cmplx(1E0_dp,0E0_dp,kind=dp)
	        end select
	    end subroutine make_poly
	    !************************************************
	    !                       max_rel_err             *
	    !************************************************
	    ! Compute maximum relative error between roots
	    ! and exact roots.
	    !************************************************
	    function max_rel_err(roots, exact_roots, deg) result(res)
	        implicit none
	        ! argument variables
	        integer, intent(in)             :: deg
	        complex(kind=dp), intent(in)    :: exact_roots(:), roots(:)
	        ! local variables
	        logical                         :: m(deg)
	        integer                         :: i, j, k
	        real(kind=dp)                   :: delta, err, res, relerr, rmax
	        ! intrinsic procedures
	        intrinsic                       :: abs, max, huge

	        ! initialize
	        m = .True.
	        rmax = huge(1.0_dp)
	        res = mu
	        ! main loop
	        do i=1,deg
	            ! find corresponding exact root (yet to be matched)
	            err = rmax
	            do j=1,deg
	                if(m(j)) then
	                    delta = abs(exact_roots(j) - roots(i))
	                    if(delta < err) then
	                        err = delta
	                        k = j
	                    end if
	                end if
	            end do
	            ! mark corresponding root as matched
	            m(k) = .False.
	            ! calculate relative error on this root and update max
	            ! zero roots give unhelpful relative errors, so revert to absolute error
	            if(abs(roots(i))<mu .or. abs(exact_roots(k))<mu) then
	                relerr = err
	            else
	                relerr = err/abs(exact_roots(k))
	            end if
	            res = max(relerr, res)
	        end do
	        res = max(res,mu)
	    end function max_rel_err
end program forwerr_multi

max_stars_repo_name: trcameron/FPML-Comp

max_issues_repo_name: trcameron/FPML-Comp

max_forks_repo_name: trcameron/FPML-Comp


text:
*> \brief \b SPOTRF2
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       RECURSIVE SUBROUTINE SPOTRF2( UPLO, N, A, LDA, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            INFO, LDA, N
*       ..
*       .. Array Arguments ..
*       REAL               A( LDA, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SPOTRF2 computes the Cholesky factorization of a real symmetric
*> positive definite matrix A using the recursive algorithm.
*>
*> The factorization has the form
*>    A = U**T * U,  if UPLO = 'U', or
*>    A = L  * L**T,  if UPLO = 'L',
*> where U is an upper triangular matrix and L is lower triangular.
*>
*> This is the recursive version of the algorithm. It divides
*> the matrix into four submatrices:
*>
*>        [  A11 | A12  ]  where A11 is n1 by n1 and A22 is n2 by n2
*>    A = [ -----|----- ]  with n1 = n/2
*>        [  A21 | A22  ]       n2 = n-n1
*>
*> The subroutine calls itself to factor A11. Update and scale A21
*> or A12, update A22 then call itself to factor A22.
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  Upper triangle of A is stored;
*>          = 'L':  Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is REAL array, dimension (LDA,N)
*>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
*>          N-by-N upper triangular part of A contains the upper
*>          triangular part of the matrix A, and the strictly lower
*>          triangular part of A is not referenced.  If UPLO = 'L', the
*>          leading N-by-N lower triangular part of A contains the lower
*>          triangular part of the matrix A, and the strictly upper
*>          triangular part of A is not referenced.
*>
*>          On exit, if INFO = 0, the factor U or L from the Cholesky
*>          factorization A = U**T*U or A = L*L**T.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*>          > 0:  if INFO = i, the leading minor of order i is not
*>                positive definite, and the factorization could not be
*>                completed.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2015
*
*> \ingroup realPOcomputational
*
*  =====================================================================
      RECURSIVE SUBROUTINE SPOTRF2( UPLO, N, A, LDA, INFO )
*
*  -- LAPACK computational routine (version 3.6.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2015
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, N
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO=0.0E+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            UPPER
      INTEGER            N1, N2, IINFO
*     ..
*     .. External Functions ..
      LOGICAL            LSAME, SISNAN
      EXTERNAL           LSAME, SISNAN
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMM, SSYRK, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      INFO = 0
      UPPER = LSAME( UPLO, 'U' )
      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SPOTRF2', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
*
*     N=1 case
*
      IF( N.EQ.1 ) THEN
*
*        Test for non-positive-definiteness
*
         IF( A( 1, 1 ).LE.ZERO.OR.SISNAN( A( 1, 1 ) ) ) THEN
            INFO = 1
            RETURN
         END IF
*
*        Factor
*
         A( 1, 1 ) = SQRT( A( 1, 1 ) )
*
*     Use recursive code
*
      ELSE
         N1 = N/2
         N2 = N-N1
*
*        Factor A11
*
         CALL SPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO )
         IF ( IINFO.NE.0 ) THEN
            INFO = IINFO
            RETURN
         END IF
*
*        Compute the Cholesky factorization A = U**T*U
*
         IF( UPPER ) THEN
*
*           Update and scale A12
*
            CALL STRSM( 'L', 'U', 'T', 'N', N1, N2, ONE,
     $                  A( 1, 1 ), LDA, A( 1, N1+1 ), LDA )
*
*           Update and factor A22
*
            CALL SSYRK( UPLO, 'T', N2, N1, -ONE, A( 1, N1+1 ), LDA,
     $                  ONE, A( N1+1, N1+1 ), LDA )
            CALL SPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
            IF ( IINFO.NE.0 ) THEN
               INFO = IINFO + N1
               RETURN
            END IF
*
*        Compute the Cholesky factorization A = L*L**T
*
         ELSE
*
*           Update and scale A21
*
            CALL STRSM( 'R', 'L', 'T', 'N', N2, N1, ONE,
     $                  A( 1, 1 ), LDA, A( N1+1, 1 ), LDA )
*
*           Update and factor A22
*
            CALL SSYRK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA,
     $                  ONE, A( N1+1, N1+1 ), LDA )
            CALL SPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
            IF ( IINFO.NE.0 ) THEN
               INFO = IINFO + N1
               RETURN
            END IF
         END IF
      END IF
      RETURN
*
*     End of SPOTRF2
*
      END

ext: f

max_stars_repo_name: sydha/Lapack

max_issues_repo_name: sydha/Lapack

max_forks_repo_name: sydha/Lapack


text:
The Plant Science Teaching Facility is a collective of buildings to support the learning of various plant sciences. It is part of the University Extension Center

Buildings
  Plant Reproductive Biology
  Bowley Hall
  The Core Greenhouse Complex
  Student Experimental Farm
  Plant Transformation Facility

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	PROGRAM MM_EVOLUTION

c	Computes apparent magnitude vs z through several filters
c	The user selects the filter(s) to use from the list in filters.log
c	The user enters values for the cosmological parameters

	include 'read_ised.dec'
	include 'filter.dec'
	include 'cosmo.dec'
	include 'SSP_13.dec'

c	Variables
	character file*128,name*128,ext*24,argu*1024
	integer kf(50)
	real yz(imw),xaux(24)
	data lun/26/,file/' '/,kf/50*0/
c	Cosmology:
	data h/70./,omega/0.30/,omega_lambda/0.70/

c	Check if correct filter file is in use.
	j=ifilter()
	if (j.eq.0) then
		write (6,'(2a)') char(7),'Please assign correct filter file'
		stop
	endif

c	Welcome user
	write (6,*)
	write (6,'(x,a)') 'Galaxy Spectral Evolution Library (GALAXEV)'
	write (6,'(x,a)') 'UNIX Version (C) 1995-2013 - G. Bruzual and S. Charlot - All Rights Reserved'

c	Ask for Ho, Omega, and Omega_lambda
	write (6,'(/x,a,f4.0,a,f5.3,a,$)') 'Enter cosmological parameters'
3	write (6,'(x,a,f4.0,a,f5.3,a,f5.3,a,$)')  'Enter Ho [',h,'], Omega [',omega,'], Omega_lambda [',omega_lambda,'] = '
	call nread(xaux,nc,*9,*10)
	if (nc.gt.0.) then
		h=xaux(1)
		omega=xaux(2)
		omega_lambda=xaux(3)
	endif
c	Omega is now entered bu the user
c	omega=1.-omega_lambda

c	Obtain cosmological constant and q
	clambda=cosmol_c(h,omega,omega_lambda,q)

c	Age of the universe and maximum age for galaxies (rounded off)
	tu=t(h,q,0.,clambda)
	tg=int(tu)*1.E9

c	Ask for galaxy age
	write (6,*)
	write (6,'(x,a,f6.2,a  )') 'Age of this universe  = tu =',tu,' Gyr'
2	write (6,'(x,a,f6.2,a,$)') 'Enter age of galaxy today = tg [',tg*1.e-9,' Gyr] = '
	call qread(ttg,nc,*8,*10)
	if (nc.gt.0.) tg=ttg*1.E9
	ttg=tg*1.E-9
	zf=zx(ttg,h,q,clambda)

c	Report cosmology
	write (6,'(/x,a)') 'Cosmological model in use:'
	write (6,100) h,omega,omega_lambda,q,clambda
100	format (' Ho = ',f4.0,3x,'Omega =',f5.2,3x,'Omega_lambda =',f5.2,3x,'qo = ',
     *            f7.4,3x,'Lambda =',1pE10.3)
	write (6,101) tu,ttg,zf
101	format (' tu = ',f5.2,' Gyr',3x,'tg = ',f5.2,' Gyr', 3x,' zf = ',f6.2)

c	Ask for *.ised file name
	write (6,*)
1	write (6,'(x,3a,$)') 'BC_GALAXEV model sed in file [',file(1:largo(file)),'] = '
	read (5,'(a)',end=10) name
	call chaext(name,'ised',nn)
	file=name
c	Read *.ised file
	call read_ised(file,tg,kl,kerr)
	if (kerr.ne.0) goto 1

c	Read filter file
	call filter0
c	Read A0Vsed
	call readA0V

c	Ask for filter numbers
4	write (6,'(/x,a)')   'Enter up to 50 filters separated by commas.'
	write (6,'( x,a,i3,a,$)') 'Compute magnitude through filter numbers (1:',nf,') = '
	read (5,'(a)',end=10) argu
	ngu = nargu(argu)
	read (argu,*,err=4) (kf(j),j=1,ngu)
	do i=1,ngu
	if (kf(i).le.0.or.kf(i).gt.nf) then
		write (6,'(x,a,i5,a)') 'Error in filter number',kf(i),'. Please try again'
		goto 4
	endif
	enddo
	write (6,'(/x,a)') 'Selected filters: '
	do i=1,ngu
	write (6,'(x,i5,x,a)') kf(i),fid(kf(i))
	enddo

c	Build output file name and write header into it
	ext='multi_mag_vega'
	call chaext(name,ext,nn)
	open (lun,file=name,form='formatted',status='unknown')
	call file_header(lun,file,1)
	write (6,'(/x,2a)') 'Output written to file(s): ',name(1:largo(name))
	ext='multi_mag_AB'
	call chaext(name,ext,nn)
	open (lun+1,file=name,form='formatted',status='unknown')
	call file_header(lun+1,file,1)
	write (6,'( x,2a)') '                           ',name(1:largo(name))

c	Compute magnitudes
	write (6,'(/x,a)') 'Computing magnitudes/colors...'
	icall=0
	do kz=1,120000
	z=zk(kz,zf)
	call evol_ised(kl,z,h,q,clambda,yz,icall,last)
	if (last.ne.0) goto 10
       	call ke_nfilt_correction(kf,w,yz,iw,z,h,q,lun,ngu)
	if (z.eq.zf) goto 10
	enddo
	close (lun)
	goto 10
!7	write (6,'(x,a)') 'Error in filter number. Please try again'
!	goto 4
8	write (6,'(x,a)') 'Error in galaxy age. Please try again'
	goto 2
9	write (6,'(x,a)') 'Error in parameters. Please try again'
	goto 3
10	end

max_stars_repo_name: jchavesmontero/galaxpy

max_issues_repo_name: jchavesmontero/galaxpy

max_forks_repo_name: jchavesmontero/galaxpy


text:
c
c
c
c    comp_R1111 -> R, dRdA, dRdf, dRdDg
c    comp_R1000 -> R
c    comp_R0111 ->    dRdA, dRdf, dRdDg
c    comp_R1100 -> R, dRdA
c    comp_R0100 ->    dRdA
c
c
c==============================================================================

      subroutine comp_R1111(R,dRdA,dRdf,dRdDg,f,A,invA,diA,An,invAn,Dg,
     *                      nip)
c
c     compute R, dRdA, dRdf, dRdDg
c
c     nip = 0 -> simple backward Euler
c         > 0 -> backward Euler exponential map with nip series terms
c
      implicit none
      integer          nip, i
      double precision R(6), dRdf(6,6), dRdA(6,6), dRdDg(6),
     *                 f(6), A(6), invA(6), diA(6,6), An(6), invAn(6),
     *                 Dg, dRdiA(6,6), hs(6), hs2(6), hu(3,3), hu2(3,3),
     *                 h4(6,6), h42(6,6), h4df(6,6), h4diA(6,6),
     *                 fact, fact2

      if (nip.eq.0) then

c       simple backward Euler
c    ------------------------------
        fact = - Dg
        call set_s_p(R,A)        ! R
        call set_s_am(R,An)      !
        call set_s_af(R,f,fact)  !
        call unit_4_f(dRdf,fact) ! dRdf
        call unit_4_p(dRdA)      ! dRdA
        call set_s_m(dRdDg,f)    ! dRdDg

      else

c       backward Euler exponential map
c    -------------------------------------
c       for i=1

        call set_s_p(R,A)
        call mult_u_ss_p(hu,invAn,A)
        call mult_s_su_am(R,A,hu)
        call set_s_p(hs,f)
        call set_s_af(R,hs,Dg)

        call unit_4_p(dRdA)
        call diff_TAT_am(dRdA,hu)

        call unit_4_f(dRdf,Dg)

        call set_s_p(dRdDg,f)

c       for i=2,3,..nip
        if (nip.gt.1) then

          call mult_u_ss_p(hu,invA,f)
          call diff_TU_p(h4,hu)
          call pzero(dRdiA,36)
          call set_4_p(h4df,h4)

c         for i=2

          fact  = Dg * Dg * 0.5d0

c         dR/df
          call mult_u_ss_p(hu2,hs,invA)
          call diff_UT_ap(h4df,hu2)
          call set_4_af(dRdf,h4df,fact)

c         dR/dinvA
          call diff_ATB_p(h4diA,hs,f)
          call set_4_af(dRdiA,h4diA,fact)

c         R
          call mult_s_su_p(hs2,hs,hu)
          call set_s_p(hs,hs2)
          call set_s_af(R,hs,fact)

c         dR/dDg
          call set_s_af(dRdDg,hs,Dg)
          fact2 = fact

c         for i=3,4...nip
          do i=3, nip

            fact  = fact * Dg / real(i)

c           dR/df
            call mult_4_44_p(h42,h4,h4df)
            call set_4_p(h4df,h42)
            call mult_u_ss_p(hu2,hs,invA)
            call diff_UT_ap(h4df,hu2)
            call set_4_af(dRdf,h4df,fact)

c           dR/dinvA
            call mult_4_44_p(h42,h4,h4diA)
            call set_4_p(h4diA,h42)
            call diff_ATB_ap(h4diA,hs,f)
            call set_4_af(dRdiA,h4diA,fact)

c           R
            call mult_s_su_p(hs2,hs,hu)
            call set_s_p(hs,hs2)
            call set_s_af(R,hs,fact)

c           dR/dDg
            call set_s_af(dRdDg,hs,fact2)
            fact2 = fact2 * Dg / real(i)

          enddo

c         dR/dA = dR/dinvA : dinvA/dA
          call mult_4_44_ap(dRdA,dRdiA,diA)

        endif

      endif

      return

      end

c=============================================================================

      subroutine comp_R1000(R,f,A,invA,An,invAn,Dg,nip)
c
c     compute R
c
c     nip = 0 -> simple backward Euler
c         > 0 -> backward Euler exponential map with nip series terms
c
      implicit none
      integer          nip, i
      double precision R(6),
     *                 f(6), A(6), invA(6), An(6), invAn(6), Dg,
     *                 hs(6), hs2(6), hu(3,3), fact

      if (nip.eq.0) then

c       simple backward Euler
c    ------------------------------
        fact = - Dg
        call set_s_p(R,A)        ! R
        call set_s_am(R,An)      !
        call set_s_af(R,f,fact)  !

      else

c       backward Euler exponential map
c    -------------------------------------
c       for i=1

        call set_s_p(R,A)
        call mult_u_ss_p(hu,invAn,A)
        call mult_s_su_am(R,A,hu)
        call set_s_p(hs,f)
        call set_s_af(R,hs,Dg)

c       for i=2,3,..nip
        if (nip.gt.1) then

          call mult_u_ss_p(hu,invA,f)

          fact  = Dg

          do i=2, nip

            fact  = fact * Dg / real(i)

c           R
            call mult_s_su_p(hs2,hs,hu)
            call set_s_p(hs,hs2)
            call set_s_af(R,hs,fact)

          enddo

        endif

      endif

      return

      end

c=============================================================================

      subroutine comp_R0111(dRdA,dRdf,dRdDg,f,A,invA,diA,An,invAn,Dg,
     *                      nip)
c
c     compute dRdA, dRdf, dRdDg
c
c     nip = 0 -> simple backward Euler
c         > 0 -> backward Euler exponential map with nip series terms
c
      implicit none
      integer          nip, i
      double precision dRdf(6,6), dRdA(6,6), dRdDg(6),
     *                 f(6), A(6), invA(6), diA(6,6), An(6), invAn(6),
     *                 Dg, dRdiA(6,6), hs(6), hs2(6), hu(3,3), hu2(3,3),
     *                 h4(6,6), h42(6,6), h4df(6,6), h4diA(6,6),
     *                 fact, fact2

      if (nip.eq.0) then

c       simple backward Euler
c    ------------------------------
        fact = - Dg
        call unit_4_f(dRdf,fact) ! dRdf
        call unit_4_p(dRdA)      ! dRdA
        call set_s_m(dRdDg,f)    ! dRdDg

      else

c       backward Euler exponential map
c    -------------------------------------
c       for i=1

        call mult_u_ss_p(hu,invAn,A)
        call set_s_p(hs,f)

        call unit_4_p(dRdA)
        call diff_TAT_am(dRdA,hu)

        call unit_4_f(dRdf,Dg)

        call set_s_p(dRdDg,f)

c       for i=2,3,..nip
        if (nip.gt.1) then

          call mult_u_ss_p(hu,invA,f)
          call diff_TU_p(h4,hu)
          call pzero(dRdiA,36)
          call set_4_p(h4df,h4)

c         for i=2

          fact  = Dg * Dg * 0.5d0

c         dR/df
          call mult_u_ss_p(hu2,hs,invA)
          call diff_UT_ap(h4df,hu2)
          call set_4_af(dRdf,h4df,fact)

c         dR/dinvA
          call diff_ATB_p(h4diA,hs,f)
          call set_4_af(dRdiA,h4diA,fact)

          call mult_s_su_p(hs2,hs,hu)
          call set_s_p(hs,hs2)

c         dR/dDg
          call set_s_af(dRdDg,hs,Dg)
          fact2 = fact

c         for i=3,4...nip
          do i=3, nip

            fact  = fact * Dg / real(i)

c           dR/df
            call mult_4_44_p(h42,h4,h4df)
            call set_4_p(h4df,h42)
            call mult_u_ss_p(hu2,hs,invA)
            call diff_UT_ap(h4df,hu2)
            call set_4_af(dRdf,h4df,fact)

c           dR/dinvA
            call mult_4_44_p(h42,h4,h4diA)
            call set_4_p(h4diA,h42)
            call diff_ATB_ap(h4diA,hs,f)
            call set_4_af(dRdiA,h4diA,fact)

            call mult_s_su_p(hs2,hs,hu)
            call set_s_p(hs,hs2)

c           dR/dDg
            call set_s_af(dRdDg,hs,fact2)
            fact2 = fact2 * Dg / real(i)

          enddo

c         dR/dA = dR/dinvA : dinvA/dA
          call mult_4_44_ap(dRdA,dRdiA,diA)

        endif

      endif

      return

      end

c=============================================================================

      subroutine comp_R1100(R,dRdA,f,A,invA,diA,An,invAn,Dg,nip)
c
c     compute R, dRdA
c
c     nip = 0 -> simple backward Euler
c         > 0 -> backward Euler exponential map with nip series terms
c
      implicit none
      integer          nip, i
      double precision R(6), dRdA(6,6),
     *                 f(6), A(6), invA(6), diA(6,6), An(6), invAn(6),
     *                 Dg, dRdiA(6,6), hs(6), hs2(6), hu(3,3), hu2(3,3),
     *                 h4(6,6), h42(6,6), h4diA(6,6),
     *                 fact

      if (nip.eq.0) then

c       simple backward Euler
c    ------------------------------
        fact = - Dg
        call set_s_p(R,A)        ! R
        call set_s_am(R,An)      !
        call set_s_af(R,f,fact)  !
        call unit_4_p(dRdA)      ! dRdA

      else

c       backward Euler exponential map
c    -------------------------------------
c       for i=1

        call set_s_p(R,A)
        call mult_u_ss_p(hu,invAn,A)
        call mult_s_su_am(R,A,hu)
        call set_s_p(hs,f)
        call set_s_af(R,hs,Dg)

        call unit_4_p(dRdA)
        call diff_TAT_am(dRdA,hu)

c       for i=2,3,..nip
        if (nip.gt.1) then

          call mult_u_ss_p(hu,invA,f)
          call diff_TU_p(h4,hu)
          call pzero(dRdiA,36)

c         for i=2

          fact  = Dg * Dg * 0.5d0

c         dR/dinvA
          call diff_ATB_p(h4diA,hs,f)
          call set_4_af(dRdiA,h4diA,fact)

c         R
          call mult_s_su_p(hs2,hs,hu)
          call set_s_p(hs,hs2)
          call set_s_af(R,hs,fact)

c         for i=3,4...nip
          do i=3, nip

            fact  = fact * Dg / real(i)

c           dR/dinvA
            call mult_4_44_p(h42,h4,h4diA)
            call set_4_p(h4diA,h42)
            call diff_ATB_ap(h4diA,hs,f)
            call set_4_af(dRdiA,h4diA,fact)

c           R
            call mult_s_su_p(hs2,hs,hu)
            call set_s_p(hs,hs2)
            call set_s_af(R,hs,fact)

          enddo

c         dR/dA = dR/dinvA : dinvA/dA
          call mult_4_44_ap(dRdA,dRdiA,diA)

        endif

      endif

      return

      end

c=============================================================================

      subroutine comp_R0100(dRdA,f,A,invA,diA,An,invAn,Dg,nip)
c
c     compute R, dRdA
c
c     nip = 0 -> simple backward Euler
c         > 0 -> backward Euler exponential map with nip series terms
c
      implicit none
      integer          nip, i
      double precision dRdA(6,6),
     *                 f(6), A(6), invA(6), diA(6,6), An(6), invAn(6),
     *                 Dg, dRdiA(6,6), hs(6), hs2(6), hu(3,3), hu2(3,3),
     *                 h4(6,6), h42(6,6), h4diA(6,6),
     *                 fact

      if (nip.eq.0) then

c       simple backward Euler
c    ------------------------------
        call unit_4_p(dRdA)      ! dRdA

      else

c       backward Euler exponential map
c    -------------------------------------
c       for i=1

        call mult_u_ss_p(hu,invAn,A)
        call set_s_p(hs,f)

        call unit_4_p(dRdA)
        call diff_TAT_am(dRdA,hu)

c       for i=2,3,..nip
        if (nip.gt.1) then

          call mult_u_ss_p(hu,invA,f)
          call diff_TU_p(h4,hu)
          call pzero(dRdiA,36)

c         for i=2

          fact  = Dg * Dg * 0.5d0

c         dR/dinvA
          call diff_ATB_p(h4diA,hs,f)
          call set_4_af(dRdiA,h4diA,fact)

          call mult_s_su_p(hs2,hs,hu)
          call set_s_p(hs,hs2)

c         for i=3,4...nip
          do i=3, nip

            fact  = fact * Dg / real(i)

c           dR/dinvA
            call mult_4_44_p(h42,h4,h4diA)
            call set_4_p(h4diA,h42)
            call diff_ATB_ap(h4diA,hs,f)
            call set_4_af(dRdiA,h4diA,fact)

            call mult_s_su_p(hs2,hs,hu)
            call set_s_p(hs,hs2)

          enddo

c         dR/dA = dR/dinvA : dinvA/dA
          call mult_4_44_ap(dRdA,dRdiA,diA)

        endif

      endif

      return

      end

c=============================================================================

ext: for

max_stars_repo_name: chennachaos/stabfem

max_issues_repo_name: chennachaos/stabfem

max_forks_repo_name: chennachaos/mpap


! NetCDF Diag Resource file library

module nc_diag_res
    ! Library to read a resource file and check if a variable is
    ! enabled within the resource file.
    !
    ! This library reads a JSON resource file with the following format:
    !   {
    !     "variables" : {
    !       "some_var" : true,
    !       "more_var" : false
    !     }
    !   }
    !
    ! Based on this sample file, we can check whether a certain variable
    ! is enabled or not using this library:
    !
    !     call nc_diag_load_resource_file("resource.json")
    !     ! This will return true:
    !     if (nc_diag_load_check_variable("some_var")) then
    !         print *, "This variable exists!"
    !         ! Do some other things here
    !     end if
    !     ! This will return false:
    !     if (nc_diag_load_check_variable("more_var")) then
    !         print *, "This variable exists!"
    !         ! Do some other things here
    !     end if
    !     ! Note that we can specify non-existent variables - these
    !     ! will also return false.
    !     if (nc_diag_load_check_variable("hmmm_var")) then
    !         print *, "This variable exists!"
    !         ! Do some other things here
    !     end if
    !     call nc_diag_close_resource_file

    use ncdres_climsg, only: ncdres_error
    use nc_diag_fson, only: ncdf_value, ncdf_parse, &
        ncdf_get, ncdf_destroy

    implicit none

    type(ncdf_value), pointer :: nc_diag_json => null()

    contains
        ! Opens a given resource file for reading.
        !
        ! Given the resource file name, open the file and set everything
        ! up for reading the file. This includes any internal memory
        ! allocation required for reading the resource file.
        !
        ! In order for memory allocation to be freed, the
        ! subroutine nc_diag_close_resource_file MUST be called.
        !
        ! If a resource file is already open, this will raise an error
        ! and the program will terminate.
        !
        ! Args:
        !     filename (character(len=*)): resource file name to load.
        !
        ! Raises:
        !     Resource file already open error if there is already a
        !     resource file currently open.
        !
        subroutine nc_diag_load_resource_file(filename)
            character(len=*), intent(in) :: filename

            if (associated(nc_diag_json)) &
                call ncdres_error("Resource file already open!")

            nc_diag_json => ncdf_parse(filename)
        end subroutine nc_diag_load_resource_file

        ! Lookup a variable and check its status.
        !
        ! Given the variable name, lookup its status within the JSON
        ! resource file.
        !
        ! If the variable is present in the JSON file, and it is
        ! enabled, this will return true. Otherwise, if the variable
        ! doesn't exist in the resource file, or it is disabled,
        ! this will return false.
        !
        ! Args:
        !     var_name (character(len=*)): variable name to lookup
        !         within the resource file.
        !
        ! Returns:
        !     var_enabled (logical): whether the variable is enabled or
        !         not within the resource file.
        !
        function nc_diag_load_check_variable(var_name) result(var_enabled)
            character(len=*), intent(in)  :: var_name
            logical                       :: var_enabled

            character(len=1024)           :: var_str

            write (var_str, "(A)") "variables." // var_name

            var_enabled = .FALSE.

            call ncdf_get(nc_diag_json, trim(var_str), var_enabled)
        end function nc_diag_load_check_variable

        ! Closes the current resource file.
        !
        ! Closes a previously opened resource file. This will free any
        ! resources allocated towards the previous resource file, and
        ! allow for opening a new resource file.
        !
        ! If no file has been opened previously, or if the file is
        ! already closed, this will raise an error and the program will
        ! terminate.
        !
        ! Raises:
        !     No resource file open error will occur if there is no
        !     resource file currently open.
        !
        subroutine nc_diag_close_resource_file
            if (associated(nc_diag_json)) then
                call ncdf_destroy(nc_diag_json)
                nullify(nc_diag_json)
            else
                call ncdres_error("No resource file open!")
            end if
        end subroutine nc_diag_close_resource_file
end module nc_diag_res


text:
      real*8 function arrmxn(array,nmax)
c
c     This subroutine finds the max norm of the real*8 array "array".
c     Note that array could actually have more than one dimension. For
c     example, if the array is really "a", which has dimensions of i1
c     and i2, this subroutine could be used by calling it in the
c     following manner: call arrmxn(a,(i1*i2)). Use this subroutine
c     with caution if nmax is not the product of the true (declared)
c     dimensions.
c
c     To find the index of the first element corresponding to the max
c     norm, call EQLIBU/iarmxn.f instead. Then just use the index to
c     get the max norm if this is also desired.
c
c     This subroutine is called by:
c
c       Any
c
c-----------------------------------------------------------------------
c
c     Input:
c
c       array  = the array
c       nmax   = the dimension or pseudo-dimension of the array
c
c     Output:
c
c       arrmxn = the max norm of the array
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nmax
c
      real*8 array(nmax)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,ileft
c
      real*8 aax,ax
c
c-----------------------------------------------------------------------
c
      arrmxn = 0.
      if (nmax .le. 0) go to 999
c
c     Note the use of a local variable (ax) within the loop.
c     Note also that the loop is unrolled.
c
      ax = 0.
      ileft = (nmax/8)*8
c
      do i = 1,ileft,8
        aax = abs(array(i))
        if (aax .gt. ax) ax = aax
        aax = abs(array(i + 1))
        if (aax .gt. ax) ax = aax
        aax = abs(array(i + 2))
        if (aax .gt. ax) ax = aax
        aax = abs(array(i + 3))
        if (aax .gt. ax) ax = aax
        aax = abs(array(i + 4))
        if (aax .gt. ax) ax = aax
        aax = abs(array(i + 5))
        if (aax .gt. ax) ax = aax
        aax = abs(array(i + 6))
        if (aax .gt. ax) ax = aax
        aax = abs(array(i + 7))
        if (aax .gt. ax) ax = aax
      enddo
c
      do i = ileft + 1,nmax
        aax = abs(array(i))
        if (aax .gt. ax) ax = aax
      enddo
c
      arrmxn = ax
c
  999 continue
      end


text:
The Wednesday Night Time Trials are a series of 10mile bicycling cycling time trials sponsored by Wheelworks and the Mad Cows Racing Team.  Held everyother Wednesday evening from Spring until Fall, the time trials are an individual race against the clock to see how quickly each participant can complete the 10mile course.  The time trials are free and open to everyone, though you must fill out a waiver to race and try really really hard not to get hit by a car.  No preregistration necessary; just be on time!

 Rules
  The course is open to traffic, so it is essential to stay to the right shoulder and yield and be courteous to vehicles.
  While waiting at the start, please do not obstruct the roadway to through traffic.
  Bicycle Helmets Bike helmets are required.
  Participants must sign a USCF waiver prior to their first time trial.
  Participants must be on time.
  Participants are sent out every 30seconds, ordered on a firstcome, firstserve basis.
  No drafting.
  Shout your number as you cross the finish line.

 Course Info

 Putah Creek Course
The Address(Putah Creek Course,38.530345,121.854687) is a flat, tenmile outandback course along Putah Creek Road between Davis and Winters.  The start/finish line is roughly threehundred yards west of Stevenson Bridge Road.  Participants will start from this point and head west towards Winters.  At the fivemile point, participants will turn around and finish back where they started.  Both the start/finish and turnaround points are marked on the road surface with permanent orange spray paint.  Warnings are also painted approaching the turnaround point, and someone is stationed at the turnaround during the time trials to make sure that participants do not miss it.

If you drive, it is recommended that you park along Stevenson Bridge Road to minimize the disruption of local traffic as much as possible.

 Bulkley Course
A flat 10 mile outandback course on Bulkley Road south of the Davis city limits. The start/finish line is on Bulkley Road, near Tremont Road. The course continues five miles south on Bulkley Road to the turnaround and returns five miles north on Bulkley back to the start/finish line.

It is recommended that you park your cars on the shoulders of Tremont Road.

 2009 Schedule
Every other Wednesday at 6 pm (on the Putah Creek course). The first time trial was on Wednesday, May 25, 6:00 PM on the Putah Creek Course.

More infos: http://stuzio.com/StuzioArticle.aspx?articleID1124

 Course Records  Putah Creek Course
Full results can be found on the Mad Cows website: http://www.madcowsracing.org/training/top250tt.php and http://www.madcowsracing.org/training/07ttresults.php .

  Mens Record  19:31  30.74 mph (http://www.paulmach.com Paul Mach, April 22, 2009)
  Womens Record (?)  22:57  26.14 mph (Molly Van Houweling, March 25, 2009)
  Team Time Trial Record  19:55  30.13 mph (Stephen Dey, Shawn Rosenthal, Tyler Dibble, Jim Wingert (didnt finish), Zach Kimsey (didnt finish); April 25, 2007)


text:
great(q1,x1).
great(c1,w1).
great(r1,v1).
great(g1,r1).
great(x1,aa1).
great(j1,v1).
great(h1,b1).
great(b1,o1).
great(q1,t1).
great(e1,r1).
great(hh1,o1).
great(z1,u1).
great(jj1,t1).
great(jj1,m1).
great(jj1,s1).
great(bb1,m1).
great(i1,aa1).
great(i1,ll1).
great(kk1,x1).
great(ff1,b1).
great(ii1,d1).
great(c1,i1).
great(k1,t1).
great(y1,u1).
great(r1,m1).
great(hh1,w1).
great(i1,o1).
great(k1,b1).
great(q1,i1).
great(l1,t1).
great(l1,u1).
great(e1,dd1).
great(g1,ff1).
great(k1,w1).
great(y1,m1).
great(d1,kk1).
great(hh1,v1).
great(ee1,bb1).
great(f1,v1).
great(kk1,t1).
great(dd1,b1).
great(c1,ee1).
great(kk1,aa1).
great(c1,t1).
great(l1,ee1).
great(ff1,x1).
great(jj1,ll1).
great(z1,t1).
great(e1,aa1).
great(v1,ee1).
great(d1,a1).
great(y1,aa1).
great(c1,y1).
great(h1,dd1).
great(k1,l1).
great(r1,o1).
great(hh1,e1).
great(f1,n1).
great(v1,aa1).
great(cc1,m1).
great(f1,ii1).
great(ff1,t1).
great(h1,x1).
great(c1,h1).
great(jj1,ee1).
great(c1,cc1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!residual for Laplace eqn!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine  block_Laplace(xloc,I_fluid,p_inter,w_inter,ien,lp_el,count_el,source)
  use global_constants
  use run_variables
  use fluid_variables
  implicit none

  integer ien(nen,ne)
  real(8) xloc(nsd,nn)
  real(8) I_fluid(nn)
  real(8) p_inter(nn)
  real(8) w_inter(nn)
!----------------------------
  integer count_el
  integer lp_el(count_el)
  real(8) source(nsd,nn)
!----------------------------
  real(8) x(nsd,nen),d(nen)
  real(8) eft0,det,effd,effm,effc
  real(8) sh(0:nsd,nen),ph(0:nsd,nen)
  real(8) xr(nsd,nsd),cf(nsd,nsd),sx(nsd,nsd)

  real(8) dr(nsd)

  integer inl,ie,isd,node,iq
!---------------------------------
  integer je
  real(8) local_source(nsd,nen)
  real(8) dg_l(nsd)
  real(8) divg
!==============================================================
  do je=1,count_el     !loop over elements
	ie=lp_el(je)
!   do ie=1,ne
     do inl=1,nen
	x(1:nsd,inl) = xloc(1:nsd,ien(inl,ie))
	d(inl) = I_fluid(ien(inl,ie))
	local_source(1:nsd,inl)=source(1:nsd,ien(inl,ie))
     end do

     do iq=1,nquad  ! loop over the quadrature points in each element
!!!!!!!!!!!!!!!!!!!!calculate shape function!!!!!!!!!!!!!!!!!!!!!!!!!!
	if (nsd==2) then
	   if(nen.eq.3) then
		include 'sh2d3n.h'
	   elseif(nen.eq.4) then
		include 'sh2d4n.h'
	   end if
	elseif(nsd==3) then
	   if(nen.eq.4) then
		include 'sh3d4n.h'
	   elseif(nen.eq.8) then
		include 'sh3d8n.h'
	   end if
	end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     eft0 = abs(det)*wq(iq)   ! calculate weight at each quad point

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     dr(:) = 0.0
     dg_l(:)=0.0

     do inl=1,nen
	dr(1:nsd) = dr(1:nsd)+sh(1:nsd,inl)*d(inl)
	dg_l(1:nsd) = dg_l(1:nsd)+sh(1:nsd,inl)*local_source(1:nsd,inl) ! dG/dx
     end do
	divg=sum(dg_l) ! div(G)

     ph(0:nsd,1:nen) = sh(0:nsd,1:nen)*eft0

     do inl=1,nen
	node=ien(inl,ie)
	do isd=1,nsd
	   p_inter(node)=p_inter(node)-ph(isd,inl)*dr(isd)+ph(0,inl)*divg
	end do! residual
	do isd=1,nsd
	   w_inter(node)=w_inter(node)+ph(isd,inl)*sh(isd,inl)
	end do! diag
     end do

     end do ! end of quad loop

  end do ! end of ele loop
!write(*,*)'w_inter=',w_inter(:)
end subroutine block_Laplace

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


text:
!****p* Program/H5test_kind
!
! NAME
!  Executable: H5test_kind
!
! FILE
!  fortran/src/H5test_kind.f90
!
! PURPOSE
!  This stand alone program is used at build time to generate the program
!  H5fortran_detect.f90. It cycles through all the available KIND parameters for
!  integers and reals. The appropriate program and subroutines are then generated
!  depending on which of the KIND values are found.
!
! NOTES
!  This program is depreciated in favor of H5test_kind_SIZEOF.f90 and is only
!  used when the Fortran intrinsic function SIZEOF is not available. It generates
!  code that does not make use of SIZEOF in H5fortran_detect.f90 which is less
!  portable in comparison to using SIZEOF.
!
!  The availability of SIZEOF is checked at configure time and the TRUE/FALSE
!  condition is set in the configure variable "FORTRAN_HAVE_SIZEOF".
!
! COPYRIGHT
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  Copyright by The HDF Group.                                               *
!  Copyright by the Board of Trustees of the University of Illinois.         *
!  All rights reserved.                                                      *
!  *
!  This file is part of HDF5.  The full HDF5 copyright notice, including     *
!  terms governing use, modification, and redistribution, is contained in    *
!  the files COPYING and Copyright.html.  COPYING can be found at the root   *
!  of the source code distribution tree; Copyright.html can be found at the  *
!  root level of an installed copy of the electronic HDF5 document set and   *
!  is linked from the top-level documents page.  It can also be found at     *
!  http://hdfgroup.org/HDF5/doc/Copyright.html.  If you do not have          *
!  access to either file, you may request a copy from help@hdfgroup.org.     *
!  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
! AUTHOR
!  Elena Pourma
!
!*****

PROGRAM test_kind
  IMPLICIT NONE
  INTEGER :: i, j, ii, ir, last, ikind_numbers(10), rkind_numbers(10)
  INTEGER :: ji, jr, jd
  last = -1
  ii = 0
  j = SELECTED_INT_KIND(18)
  DO i = 1,100
     j = SELECTED_INT_KIND(i)
     IF(j .NE. last) THEN
        IF(last .NE. -1) THEN
           ii = ii + 1
           ikind_numbers(ii) = last
        ENDIF
        last = j
        IF(j .EQ. -1) EXIT
     ENDIF
  ENDDO

  last = -1
  ir = 0
  DO i = 1,100
     j = SELECTED_REAL_KIND(i)
     IF(j .NE. last) THEN
        IF(last .NE. -1) THEN
           ir = ir + 1
           rkind_numbers(ir) = last
        ENDIF
        last = j
        IF(j .EQ. -1) EXIT
     ENDIF
  ENDDO

!  Generate program information:

WRITE(*,'(40(A,/))') &
'!****h* ROBODoc/H5fortran_detect.f90',&
'!',&
'! NAME',&
'!  H5fortran_detect',&
'! ',&
'! PURPOSE',&
'!  This stand alone program is used at build time to generate the header file',&
'!  H5fort_type_defines.h. The source code itself was automatically generated by',&
'!  the program H5test_kind.f90',&
'!',&
'! NOTES',&
'!  This source code does not make use of the Fortran intrinsic function SIZEOF because',&
'!  the availability of the intrinsic function was determined to be not available at',&
'!  configure time',&
'!',&
'! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *',&
'!   Copyright by The HDF Group.                                               *',&
'!   Copyright by the Board of Trustees of the University of Illinois.         *',&
'!   All rights reserved.                                                      *',&
'!                                                                             *',&
'!   This file is part of HDF5.  The full HDF5 copyright notice, including     *',&
'!   terms governing use, modification, and redistribution, is contained in    *',&
'!   the files COPYING and Copyright.html.  COPYING can be found at the root   *',&
'!   of the source code distribution tree; Copyright.html can be found at the  *',&
'!   root level of an installed copy of the electronic HDF5 document set and   *',&
'!   is linked from the top-level documents page.  It can also be found at     *',&
'!   http://hdfgroup.org/HDF5/doc/Copyright.html.  If you do not have          *',&
'!   access to either file, you may request a copy from help@hdfgroup.org.     *',&
'! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *',&
'!',&
'! AUTHOR',&
'!  H5test_kind.f90',&
'!',&
'!*****'

!  Generate a program

  WRITE(*,*) "PROGRAM int_kind"
  WRITE(*,*) "WRITE(*,*) "" /*generating header file*/ """
  ji = 0
  WRITE(*, "("" CALL i"", i2.2,""()"")") ji
  jr = 0
  WRITE(*, "("" CALL r"", i2.2,""()"")") jr
  jd = 0
  WRITE(*, "("" CALL d"", i2.2,""()"")") jd
  DO i = 1, ii
     j = ikind_numbers(i)
     WRITE(*, "("" CALL i"", i2.2,""()"")") j
  ENDDO
  DO i = 1, ir
     j = rkind_numbers(i)
     WRITE(*, "("" CALL r"", i2.2,""()"")") j
  ENDDO
  WRITE(*,*) "END PROGRAM int_kind"
  j = 0
  ji = KIND(1)
  WRITE(*, "("" SUBROUTINE i"", i2.2,""()"")") j
  WRITE(*,*)"   IMPLICIT NONE"
  WRITE(*,*)"   INTEGER :: a = 0"
  WRITE(*,*)"   INTEGER :: a_size"
  WRITE(*,*)"   CHARACTER(LEN=2) :: jchr2"
  WRITE(*,*)"   a_size = BIT_SIZE(a)"
  WRITE(*,*)"   IF (a_size .EQ. 8) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",ji
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_NATIVE_1_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   endif"
  WRITE(*,*)"   IF (a_size .EQ. 16) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",ji
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_NATIVE_2_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   endif"
  WRITE(*,*)"   IF (a_size .EQ. 32) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",ji
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_NATIVE_4_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   IF (a_size .EQ. 64) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",ji
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_NATIVE_8_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   IF (a_size .EQ. 128) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",ji
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_NATIVE_16_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   RETURN"
  WRITE(*,*)"END SUBROUTINE"
  jr = KIND(1.0)
  WRITE(*, "("" SUBROUTINE r"", i2.2,""()"")") j
  WRITE(*,*)"   IMPLICIT NONE"
  WRITE(*,*)"   REAL :: b(32)"
  WRITE(*,*)"   INTEGER :: a(1)"
  WRITE(*,*)"   INTEGER :: a_size"
  WRITE(*,*)"   INTEGER :: real_size"
  WRITE(*,*)"   CHARACTER(LEN=2) :: jchr2"
  WRITE(*,*)"   a_size = BIT_SIZE(a(1)) ! Size in bits for integer"
  WRITE(*,*)"   real_size = (SIZE(TRANSFER(b,a))*a_size)/SIZE(b)"
  WRITE(*,*)"   IF (real_size .EQ. 32) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",jr
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_REAL_NATIVE_4_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   IF (real_size .EQ. 64) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",jr
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_REAL_NATIVE_8_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   IF (real_size .EQ. 128) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",jr
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_REAL_NATIVE_16_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   RETURN"
  WRITE(*,*)"END SUBROUTINE"
  jd = KIND(1.d0)
  WRITE(*, "("" SUBROUTINE d"", i2.2,""()"")") j
  WRITE(*,*)"   IMPLICIT NONE"
  WRITE(*,*)"   DOUBLE PRECISION :: b=0"
  WRITE(*,*)"   INTEGER :: a(8)=0"
  WRITE(*,*)"   INTEGER :: a_size"
  WRITE(*,*)"   INTEGER :: b_size"
  WRITE(*,*)"   CHARACTER(LEN=2) :: jchr2"
  WRITE(*,*)"   a_size = BIT_SIZE(a(1))"
  WRITE(*,*)"   b_size = SIZE(transfer(b,a))*a_size"
  WRITE(*,*)"   IF (b_size .EQ. 64) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",jd
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_DOUBLE_NATIVE_8_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   IF (b_size .EQ. 128) THEN"
  WRITE(*,*)"       WRITE(jchr2,'(I2)')",jd
  WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_DOUBLE_NATIVE_16_KIND "'//"//ADJUSTL(jchr2)"
  WRITE(*,*)"   ENDIF"
  WRITE(*,*)"   RETURN"
  WRITE(*,*)"END SUBROUTINE"
  DO i = 1, ii
     j = ikind_numbers(i)
     WRITE(*, "("" SUBROUTINE i"", i2.2,""()"")") j
     WRITE(*,*)"   IMPLICIT NONE"
     WRITE(*,*)"   INTEGER(",j,") :: a = 0"
     WRITE(*,*)"   INTEGER :: a_size"
     WRITE(*,*)"   CHARACTER(LEN=2) :: jchr2"
     WRITE(*,*)"   a_size = BIT_SIZE(a)"
     WRITE(*,*)"   IF (a_size .EQ. 8) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_INTEGER_1_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   ENDIF"
     WRITE(*,*)"   IF (a_size .EQ. 16) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_INTEGER_2_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   ENDIF"
     WRITE(*,*)"   IF (a_size .EQ. 32) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_INTEGER_4_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   ENDIF"
     WRITE(*,*)"   IF (a_size .EQ. 64) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_INTEGER_8_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   ENDIF"
     WRITE(*,*)"   IF (a_size .EQ. 128) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_INTEGER_16_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   ENDIF"
     WRITE(*,*)"   RETURN"
     WRITE(*,*)"   END SUBROUTINE"
  ENDDO
  DO i = 1, ir
     j = rkind_numbers(i)
     WRITE(*, "("" SUBROUTINE r"", i2.2,""()"")") j
     WRITE(*,*)"   IMPLICIT NONE"
     WRITE(*,*)"   REAL(KIND=",j,") :: b(32)"
     WRITE(*,*)"   INTEGER :: a(1)"
     WRITE(*,*)"   INTEGER :: a_size"
     WRITE(*,*)"   INTEGER :: real_size"
     WRITE(*,*)"   CHARACTER(LEN=2) :: jchr2"
     WRITE(*,*)"   a_size = BIT_SIZE(a(1)) ! Size in bits for integer"
     WRITE(*,*)"   real_size = (SIZE(TRANSFER(b,a))*a_size)/SIZE(b)"
     WRITE(*,*)"   IF (real_size .EQ. 32) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_REAL_4_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   ENDIF"
     WRITE(*,*)"   IF (real_size .EQ. 64) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,*)'       WRITE(*,*) "#define H5_FORTRAN_HAS_REAL_8_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   endif"
     WRITE(*,*)"   IF (real_size .EQ. 128) THEN"
     WRITE(*,*)"       WRITE(jchr2,'(I2)')",j
     WRITE(*,'(A)')'       WRITE(*,*) "#define H5_FORTRAN_HAS_REAL_16_KIND "'//"//ADJUSTL(jchr2)"
     WRITE(*,*)"   ENDIF"
     WRITE(*,*)"   RETURN"
     WRITE(*,*)"   END SUBROUTINE"
  ENDDO
END PROGRAM test_kind

max_stars_repo_name: pbasting/cactus

max_issues_repo_name: pbasting/cactus

max_forks_repo_name: pbasting/cactus