*======================================================================= ! This component has been assigned Type Number 204. If that number conflicts with ! another user Type number, you will need to change it and recompile the appropriate ! dll. subroutine type204 (time, xin, out, t, dtdt, par, info, icntrl, *) !DEC$ATTRIBUTES DLLEXPORT :: TYPE204 *----------------------------------------------------------------------- * Unglazed transpired collector (UTC) subroutine. * * Reference: Kutscher, C.F., "Heat Exchanger Effectiveness and * Pressure Drop for Air Flow Through Perforated Plates With and * Without Crosswind," J. Heat Transfer, vol 116, p 391, 1994. * * Parameters: * 1. Collector area (m2) * 2. Collector height (m) * 3. Collector hole diameter (m) * 4. Collector hole pitch, distance between centers of holes (m) * 5. Collector emissivity * 6. Collector absorptivity * 7. Plenum depth (m) * 8. Emissivity of the wall behind the collector * 9. R-value of the wall behind the collector (C-m2-hr/kJ) * 10. Total UA-value of the building walls and roof (kJ/hr-C) * 11. Room air temperature (C) * 12. Ambient air temperature above which the summer bypass damper * is opened (C) * 13. Maximum auxiliary heat rate available (kJ/hr) * 14. Night bypass = 0 if bypass not automatically opened at night * = 1 if bypass automatically opened at night * * Inputs: * 1. Month of year * 2. Hour of month * 3. Radiation incident on the collector (kJ/m2) * 4. Ambient temperature (C) * 5. Sky temperature (C) * 6. Atmospheric pressure (kPa) * 7. Internal gains due to people, equipment, etc. (kJ/hr) * 8. Supply air flow rate from collector air-handling units (m3/hr) * 9. Minimum outdoor air flow rate through collector/summer bypass * damper (m3/hr) * 10. Supply air flow rate from no-collector air-handling units * (m3/hr) * 11. Outdoor air flow rate through no collector (m3/hr) * * Outputs: * 1. Plenum air temperature (C) * 2. Collector outlet air temperature (C) * 3. Mixed air temperature (C) * 4. Supply air temperature (C) * 5. Collector surface temperature (C) * 6. Energy savings rate (kJ/hr) * 7. Convection from collector (kJ/hr) * 8. Convection from wall (kJ/hr) * 9. Radiation from collector (kJ/hr) * 10. Radiation from wall (kJ/hr) * 11. Conduction through wall (kJ/hr) * 12. Reduced conduction through wall because of collector (kJ/hr) * 13. Absorbed energy rate (kJ/hr) * 14. Auxiliary heating rate (kJ/hr) * 15. Outdoor air flow rate through collector/summer bypass damper * (m3/hr) * 16. Heat exchanger effectiveness of collector (C) * 17. Solar efficiency of the collector * 18. Pressure drop across collector plate (kPa) * 19. Bypass damper position = 0.0 if open * = 1.0 if closed * 20. Heat rate supplied by a traditional heating system (kJ/hr) * 21. Additional fan power required (kJ/hr) * * For this model to correctly calculate the performance of transpired * collectors, the approach velocity (appvel) should be greater than * 72 m/hr. Otherwise, there will be convection losses from the * collector between the holes, and the collector's performance will be * reduced. Also, the collector pressure drop (pcol) should be at least * 0.025 Pa to ensure uniform flow through the collector. Otherwise, * sections of the collector will become hotter than others, and * radiation losses from the collector will increase. Again, this will * reduce the collector's performance. To achieve a sufficient pressure * drop, the porosity (por) should be about 0.005 to 0.01 for the given * approach velocities. If these approach velocity and pressure drop * conditions are not met, this subroutine will write a warning to a * file. It is important to emphasize that a collector _can_ be operated * at approach velocities and pressure drops below these values, but this * model will just over predict the performance. * * Written by David Summers, Solar Energy Lab, U. Wisconsin - Madison * for M.S. thesis, December 1995. *----------------------------------------------------------------------- implicit none * type71 variables * parameters real*8 area, ht, diam, pitch, emisc, absor, depth, emisw, rwall real*8 ua, troom, tbypass, qmax, nitebp * inputs real*8 month, hour, rad, tamb, tsky, patm, gain, flow1, minout1 real*8 flow2, out2 * values calculated directly from parameters and inputs real*8 por, surfa, udwall, qabs, tgnd, tsur, qbldg real*8 aircp, aircond, airvisc, airden, a0, a1, a2 real*8 mflow1, mflow2, mflow, gmin, beta, qtrad * outputs and other variables real*8 gamma, tout, tcol, tplen, twall, qvcol, qrcol, qvwall real*8 qrwall, qdwall, red, effhx, appvel real*8 tmix, tsup, qaux1, lowg, hig, oldg, dif, diflim, small, g real*8 zeta, pcol, plenden, f, dh, plenvel, pfric, pbuoy, pacc real*8 delp, fanpw real*8 qaux2, qaux, film, tsolair, qpot, qred, qu, qsave, soleff integer istr, i, j, jmax real warn character*10 warnfile * TRNSYS variables character*3 ycheck, ocheck real time, t, dtdt, par, s, time0, tfinal, delt double precision xin, out integer*4 info, unit, type integer np, ni, no, icntrl, iwarn, nstore, iav, iunit integer lur, luw, iform, luk parameter ( np=14, ni=11, no=21 ) dimension par(np), xin(ni), out(no), info(15) dimension ycheck(ni), ocheck(no) * UTC common block (used in utcsolve subroutine) common /utc/ area, ht, diam, pitch, emisc, depth, emisw, troom, & flow1, mflow1, por, surfa, udwall, tamb, qabs, tsur, & aircp, aircond, airvisc, airden, unit, type, & month, hour, g * TRNSYS common blocks common /sim/ time0, tfinal, delt, iwarn common /store/ nstore, iav, s(5000) common /lunits/ lur, luw, iform, luk data iunit/0/ ! Set the version information for TRNSYS IF (INFO(7).EQ.-2) THEN INFO(12) = 15 RETURN 1 ENDIF *----------------------------------------------------------------------- diflim = 1.0e-4 small = 1.0e-4 jmax = 100 warnfile = 'WARN.UTC' g = 9.8 * 3600**2 ! acceleration of gravity, m/hr2 *----------------------------------------------------------------------- * first call of simulation if ( info(7) .gt. -1 ) go to 10 info(10) = 1 * check parameters info(6) = no call typeck( 1, info, ni, np, 0 ) * setup for warnings unit = info(1) type = info(2) istr = info(10) s(istr) = -1.0 open( unit = 71, file = warnfile, status = 'unknown' ) write (71,800) type 800 format ('Type',i3,' warnings:') * set variable types data ycheck/'MN1','TD1','IR1','TE1','TE1','PR2','PW1','VF1','VF1', & 'VF1','VF1'/ data ocheck/'TE1','TE1','TE1','TE1','TE1','PW1','PW1','PW1','PW1', & 'PW1','PW1','PW1','PW1','PW1','VF1','DM1','DM1','PR2', & 'DM1','PW1','PW1'/ call rcheck( info, ycheck, ocheck ) *----------------------------------------------------------------------- * if its a different unit, set parameters 10 continue if ( info(1) .eq. iunit ) go to 20 iunit = info(1) area = par(1) ht = par(2) diam = par(3) pitch = par(4) emisc = par(5) absor = par(6) depth = par(7) emisw = par(8) rwall = par(9) ua = par(10) troom = par(11) tbypass = par(12) qmax = par(13) nitebp = par(14) por = 0.907 * ( diam / pitch )**2 surfa = ( 1 - por ) * area udwall = 1.0 / rwall *----------------------------------------------------------------------- * set inputs 20 continue month = xin(1) hour = xin(2) rad = xin(3) tamb = xin(4) tsky = xin(5) patm = xin(6) gain = xin(7) flow1 = xin(8) minout1 = xin(9) flow2 = xin(10) out2 = xin(11) qabs = rad * surfa * absor tgnd = tamb tsur = ( 0.5 * ((tsky+273.15)**4 + (tgnd+273.15)**4) )**0.25 & - 273.15 * calculate building loss, assuming no reduced wall loss qbldg = ua * (troom - tamb) - gain * calculate ambient air properties (curve-fits from EES) aircp = 1.0062 + 3.6028e-5*tamb & - 1.0885e-6*tamb*tamb + 1.3791e-8*tamb**3 aircond = 0.08659 + 2.6993e-4*tamb airvisc = 0.06201 + 1.7428e-4*tamb a0 = -5.1936e-5 + 1.2758e-2*patm a1 = -5.7037e-7 - 4.6865e-5*patm a2 = -5.6746e-9 + 1.5511e-7*patm airden = a0 + a1*tamb + a2*tamb*tamb mflow1 = flow1 * airden mflow2 = flow2 * airden mflow = mflow1 + mflow2 qtrad = (minout1 + out2) * airden * aircp * (troom - tamb) & + qbldg if ( qtrad .lt. small ) qtrad = 0.0 *----------------------------------------------------------------------- * check if summer bypass damper is open if ( tamb .gt. tbypass ) then do i = 1,no out(i) = 0.0 enddo out(2) = tamb ! tout = tamb out(3) = tamb ! tmix = tamb out(4) = tamb ! tsup = tamb out(15) = flow1 ! 100% fresh air return 1 endif *----------------------------------------------------------------------- * check air flow rates if ( minout1 .gt. flow1 .or. out2 .gt. flow2 ) then write (luw,810) unit, type 810 format (//1x,'***** ERROR *****',/1x,'UNIT ',i3, & ' TYPE ',i3,' UNGLAZED TRANSPIRED COLLECTOR',/1x, & 'CHECK AIR FLOW RATES IN DECK') CALL MYSTOP(1001) RETURN 1 endif if ( flow2 .gt. small ) then beta = out2 / flow2 else beta = 0.0 endif if ( flow1 .gt. small ) then gmin = minout1 / flow1 else qaux = qtrad do i = 1,no out(i) = 0.0 enddo out(2) = tamb ! tout = tamb out(3) = tamb ! tmix = tamb out(4) = tamb ! tsup = tamb out(14) = qaux out(20) = qtrad return 1 endif *----------------------------------------------------------------------- * at night, open bypass damper if nitebp = 1 if ( qabs .lt. small .and. nitebp .gt. 0.5 ) then gamma = gmin tout = tamb tmix = gamma*tout + (1.0 - gamma)*troom tsup = troom + qbldg / (mflow * aircp) qaux = qtrad do i = 1,no out(i) = 0.0 enddo out(2) = tout out(3) = tmix out(4) = tsup out(14) = qaux out(15) = minout1 out(20) = qtrad return 1 endif *----------------------------------------------------------------------- * find gamma such that qaux1 is minimized. in the winter, this will * usually be when gamma = gmin. if the air is too hot for this * condition, then no auxiliary heat is needed (i.e. qaux1 = 0) and * gamma > gmin. if gamma is between gmin and 1.0, then use the * bisection method to determine gamma. *----------------------------------------------------------------------- * calculate necessary supply air temperature to meet load tsup = troom + qbldg / (mflow * aircp) * solve equations for gamma = gmin gamma = gmin call utcsolve( gamma, tout, tcol, tplen, twall, & qvcol, qrcol, qvwall, qrwall, qdwall, & red, effhx, appvel ) tmix = gamma*tout + (1.0 - gamma)*troom if ( tmix .lt. tsup ) then * calculate auxiliary heat qaux1 = mflow1 * aircp * ( tsup - tmix ) else * solve equations for gamma = 1 qaux1 = 0.0 gamma = 1.0 call utcsolve( gamma, tout, tcol, tplen, twall, & qvcol, qrcol, qvwall, qrwall, qdwall, & red, effhx, appvel ) tmix = tout if ( tmix .lt. tsup ) then * use bisection method to find gamma j = 0 lowg = gmin hig = 1.0 gamma = (lowg + hig) / 2.0 100 continue j = j + 1 call utcsolve( gamma, tout, tcol, tplen, twall, & qvcol, qrcol, qvwall, qrwall, qdwall, & red, effhx, appvel ) tmix = gamma*tout + (1.0 - gamma)*troom if ( tmix .lt. tsup ) then hig = gamma else lowg = gamma endif oldg = gamma gamma = (lowg + hig) / 2.0 dif = abs( gamma - oldg ) 109 if ( dif .gt. diflim .and. j .lt. jmax ) go to 100 if ( j .ge. jmax ) then write (luw,811) unit, type, dif, nint(month), nint(hour) 811 format (//1x,'***** ERROR *****',/1x,'UNIT ',i3, & ' TYPE ',i3,' UNGLAZED TRANSPIRED COLLECTOR',/1x, & 'NO CONVERGENCE IN J LOOP, DIF = ',e7.2,/1x, & 'MONTH = ',i2,' HOUR = ',i3) CALL MYSTOP(1001) RETURN endif endif endif *----------------------------------------------------------------------- * calculate additional fan power * calculate pressure drop across collector zeta = 6.82 * ( (1-por)/por )**2 * ( red )**( -0.236 ) pcol = 0.5 * airden * appvel * appvel * zeta * 7.71605e-11 ! factor of 7.7e-11 to convert units * calculate air density in plenum plenden = a0 + a1*tplen + a2*tplen*tplen * calculate friction pressure drop through plenum f = 0.05 ! estimate dh = 4.0 * ( depth * area/ht ) / ( 2.0 * ( area/ht + depth ) ) plenvel = 0.5 * appvel * ht / depth pfric = f * (ht / dh) * plenden * plenvel**2 / 2.0 * convert units to kPa pfric = pfric / ( 3600.0**2 * 1000.0 ) * calculate buoyancy pressure term pbuoy = (airden - plenden) * g * ht * convert units to kPa pbuoy = pbuoy / ( 3600.0**2 * 1000.0 ) * calculate acceleration pressure drop pacc = plenden * (2.0*plenvel)**2 / 2.0 * convert units to kPa pacc = pacc / ( 3600.0**2 * 1000.0 ) * calculate total pressure drop and fan power delp = pcol + pfric - pbuoy + pacc fanpw = gamma * flow1 * delp *----------------------------------------------------------------------- * calculate reduced conduction through wall film = 54.0 ! (kJ/hr-m2-C) film coefficient for air ! against original wall tsolair = tamb + (emisw * rad / film) ! sol-air temp for ! absor = emis qpot = udwall * area * (troom - tsolair) ! potential conduction qred = qpot - qdwall ! reduced conduction * calculate total auxiliary heat required, subtracting reduced wall * loss since it was assumed to be zero when qbldg was calculated qaux2 = mflow2 * aircp * ( tsup - beta*tamb - (1.0-beta)*troom ) qaux = qaux1 + qaux2 - qred if ( qaux .lt. small ) qaux = 0.0 * calculate energy savings qu = qvcol + qvwall ! useful energy gained by air qsave = qtrad - qaux ! energy saved * calculate solar efficiency if ( rad .gt. small ) then soleff = qvcol / ( rad * area ) if ( soleff .gt. 1.0 ) soleff = 1.0 if ( soleff .lt. 0.0 ) soleff = 0.0 else soleff = 0.0 endif *----------------------------------------------------------------------- * write warnings to file if ( info(7) .lt. 0 ) go to 30 warn = -1.0 if ( appvel .lt. 72.0 ) then write (71,801) month, hour 801 format (/,'month = ',f3.0,' hour = ',f5.1) write (71,802) appvel 802 format ('* approach velocity = ',f5.2,' m/hr',/,' For appvel', & ' < 72 m/hr, model overpredicts performance') warn = 1.0 endif if ( pcol .lt. 0.025 ) then if ( warn .lt. 0.0 ) then write (71,801) month, hour endif write (71,803) pcol 803 format ('* pressure drop = ',f6.4,' kPa',/,' For pcol', & ' < 0.025 kPa, model overpredicts performance') warn = 1.0 endif if ( qaux .gt. qmax ) then if ( warn .lt. 0.0 ) then write (71,801) month, hour endif write (71,804) qaux, qmax 804 format ('* qaux = ',e8.2,' kJ/hr; qmax = ',e8.2,' kJ/hr',/, & ' auxiliary heater(s) too small to meet load') warn = 1.0 endif if ( tmix .gt. tsup ) then if ( warn .lt. 0.0 ) then write (71,801) month, hour endif write (71,805) tmix, tsup 805 format ('* actual tsup = ',f5.1,' C; desired tsup = ',f5.1, & ' C',/,' tsup is too hot, summer bypass damper should', & ' have been opened') tsup = tmix warn = 1.0 endif if ( warn .gt. 0.0 ) then istr = info(10) if ( s(istr) .lt. 0.0 ) then write (luw,812) unit, type, warnfile 812 format (//2x,'***** WARNING ***** UNIT',i3,' TYPE',i3/4x, & 'CHECK ', a10, & ' FOR UNGLAZED TRANSPIRED COLLECTOR WARNINGS') s(istr) = 1.0 iwarn = iwarn + 1 endif endif 30 continue *----------------------------------------------------------------------- * set outputs out(1) = tplen out(2) = tout out(3) = tmix out(4) = tsup out(5) = tcol out(6) = qsave out(7) = qvcol out(8) = qvwall out(9) = qrcol out(10) = qrwall out(11) = qdwall out(12) = qred out(13) = qabs out(14) = qaux out(15) = gamma * flow1 out(16) = effhx out(17) = soleff out(18) = pcol out(19) = 1.0 out(20) = qtrad out(21) = fanpw *----------------------------------------------------------------------- return 1 end *======================================================================= subroutine utcsolve( gamma, tout, tcol, tplen, twall, & qvcol, qrcol, qvwall, qrwall, qdwall, & red, effhx, appvel ) *----------------------------------------------------------------------- * This subroutine solves the energy balances on the collector, air, and * outside wall surface. The temperatures and heat flows are output. * * The temperatures are: * * tout = air at the outlet from the collector * tcol = collector surface * tplen = air in the plenum * twall = outside wall surface * troom = air in the room * tamb = ambient air * tsur = surroundings (sky & ground) for radiation calculation * * All temperatures in this subroutine are converted from Celsius to * Kelvin at the beginning, and then they are converted back to Celsius * at the end for the main UTC subroutine. * The heat flows are labelled 'qXsource', where X = (r, v, d) for * radiation, convection, conduction. The source of the heat flow is * defined for the usual direction of heat flow for winter operation, * from the inside of the building to the outside. So, qvwall is the * convection from the wall to the plenum (not wall to the room). * Similarly, qrcol is the radiation from the collector to the * surroundings (not collector to wall). The only exceptions are qabs, * the absorbed solar energy, and qdwall, the conduction through the * wall. The heat transfer coefficients are labelled the same way. * * Written by David Summers, Solar Energy Lab, U. Wisconsin - Madison * for M.S. thesis, December 1995. *----------------------------------------------------------------------- implicit none * TRNSYS variables integer lur, luw, iform, luk * utcsolve variables * arguments real*8 gamma, tout, tcol, tplen, twall, qvcol, qrcol, qvwall real*8 qrwall, qdwall, red, effhx, appvel * UTC common block variables real*8 area, ht, diam, pitch, emisc, depth, emisw, troom real*8 flow1, mflow1, por, surfa, udwall, tamb, qabs, tsur real*8 aircp, aircond, airvisc, airden, month, hour, g integer*4 unit, type * subroutine internal variables real*8 pi, sb, diflim, small real*8 holevel, plenvel, nud, hvcol, pr, reht, nuht, hvwall double precision a, b real*8 x, res, dif integer i, j, k, kmax, neq, one, ipiv, ifail parameter ( neq = 8 ) dimension a(neq,neq), b(neq), x(neq), ipiv(neq), res(neq) * TRNSYS common block common /lunits/ lur, luw, iform, luk * UTC common block common /utc/ area, ht, diam, pitch, emisc, depth, emisw, troom, & flow1, mflow1, por, surfa, udwall, tamb, qabs, tsur, & aircp, aircond, airvisc, airden, unit, type, & month, hour, g *----------------------------------------------------------------------- pi = 3.14159265359 sb = 2.0412e-7 ! Stefan-Boltzmann, kJ/hr-m2-K4 diflim = 1.0e-4 small = 1.0e-4 kmax = 100 one = 1 ifail = 0 * check for no flow through collector if ( gamma .lt. small ) then tout = tamb tcol = tamb tplen = tamb twall = tamb qvcol = 0.0 qrcol = 0.0 qvwall = 0.0 qrwall = 0.0 qdwall = udwall * area * ( troom - tplen ) red = 0.0 effhx = 0.0 appvel = 0.0 return endif * convert celsius to kelvin troom = troom + 273.15 tamb = tamb + 273.15 tsur = tsur + 273.15 * calculate approach, hole, and plenum velocities appvel = gamma * flow1 / area holevel = gamma * flow1 / ( area * por ) plenvel = 0.5 * appvel * ht / depth * calculate heat exchanger effectiveness for collector red = airden * holevel * diam / airvisc nud = 2.75 * ( pitch/diam )**( -1.2 ) * ( red )**( 0.43 ) hvcol = nud * aircond / diam effhx = 1 - exp( - hvcol * surfa / (gamma * mflow1 * aircp) ) * calculate heat transfer coefficient for wall to air convection pr = 0.71 reht = airden * plenvel * ht / airvisc if ( reht .gt. 500000 ) then * turbulent nuht = (0.037*reht**0.8 - 871.0) * pr**(1.0/3.0) else * laminar nuht = 0.664 * reht**0.5 * pr**(1.0/3.0) endif hvwall = nuht * aircond / ht *----------------------------------------------------------------------- * solve simultaneous equations for unknown temperatures and heat flows * [a][x] = [b] *----------------------------------------------------------------------- * initial guesses for collector and wall temperatures tcol = tamb twall = tamb * calculate [a] matrix and [b] array do i = 1,neq do j = 1,neq a(i,j) = 0.0 enddo b(i) = 0.0 enddo a(1,1) = - effhx a(1,2) = 1.0 a(2,1) = - emisc * sb * surfa & * (tcol**2 + tsur**2) * (tcol + tsur) a(2,3) = 1.0 a(3,1) = sb * area * (tcol**2 + twall**2) * (tcol + twall) & / ( 1.0/emisw + 1.0/emisc - 1.0 ) a(3,4) = 1.0 a(3,5) = - sb * area * (tcol**2 + twall**2) * (tcol + twall) & / ( 1.0/emisw + 1.0/emisc - 1.0 ) a(4,2) = hvwall * area a(4,5) = - hvwall * area a(4,6) = 1.0 a(5,2) = - gamma * mflow1 * aircp a(5,7) = 1.0 a(6,3) = 1.0 a(6,4) = - 1.0 a(6,7) = 1.0 a(7,4) = - 1.0 a(7,6) = - 1.0 a(7,8) = 1.0 a(8,5) = udwall * area a(8,8) = 1.0 b(1) = tamb * ( 1.0 - effhx ) b(2) = - emisc * sb * surfa & * (tcol**2 + tsur**2) * (tcol + tsur) * tsur b(5) = - gamma * mflow1 * aircp * tamb b(6) = qabs b(8) = udwall * area * troom * start iterations k = 0 110 continue k = k + 1 * solve for unknowns using lapack matrix solver * see http://www.netlib.org -> lapack -> lapack/double call dgesv( neq, one, a, neq, ipiv, b, neq, ifail ) if ( ifail .ne. 0 ) then write (luw,813) unit, type, ifail, nint(month), nint(hour) 813 format (//1x,'***** ERROR *****',/1x,'UNIT ',i3, & ' TYPE ',i3,' UNGLAZED TRANSPIRED COLLECTOR',/1x, & 'MATRIX SOLVER ERROR, IFAIL = ',i3,/1x, & 'MONTH = ',i2,' HOUR = ',i3) CALL MYSTOP(1001) RETURN endif do i = 1,neq x(i) = b(i) enddo * calculate [a] matrix and [b] array tcol = x(1) twall = x(5) do i = 1,neq do j = 1,neq a(i,j) = 0.0 enddo b(i) = 0.0 enddo a(1,1) = - effhx a(1,2) = 1.0 a(2,1) = - emisc * sb * surfa & * (tcol**2 + tsur**2) * (tcol + tsur) a(2,3) = 1.0 a(3,1) = sb * area * (tcol**2 + twall**2) * (tcol + twall) & / ( 1.0/emisw + 1.0/emisc - 1.0 ) a(3,4) = 1.0 a(3,5) = - sb * area * (tcol**2 + twall**2) * (tcol + twall) & / ( 1.0/emisw + 1.0/emisc - 1.0 ) a(4,2) = hvwall * area a(4,5) = - hvwall * area a(4,6) = 1.0 a(5,2) = - gamma * mflow1 * aircp a(5,7) = 1.0 a(6,3) = 1.0 a(6,4) = - 1.0 a(6,7) = 1.0 a(7,4) = - 1.0 a(7,6) = - 1.0 a(7,8) = 1.0 a(8,5) = udwall * area a(8,8) = 1.0 b(1) = tamb * ( 1.0 - effhx ) b(2) = - emisc * sb * surfa & * (tcol**2 + tsur**2) * (tcol + tsur) * tsur b(5) = - gamma * mflow1 * aircp * tamb b(6) = qabs b(8) = udwall * area * troom * calculate residuals do i = 1,neq res(i) = 0.0 do j = 1,neq res(i) = res(i) + a(i,j) * x(j) enddo res(i) = res(i) - b(i) enddo * check for convergence dif = 0.0 do i = 1,neq dif = sqrt( dif**2 + res(i)**2 ) enddo 119 if ( dif .gt. diflim .and. k .lt. kmax ) go to 110 if ( k .ge. kmax ) then write (luw,815) unit, type, dif, nint(month), nint(hour) 815 format (//1x,'***** ERROR *****',/1x,'UNIT ',i3, & ' TYPE ',i3,' UNGLAZED TRANSPIRED COLLECTOR',/1x, & 'NO CONVERGENCE IN K LOOP, DIF = ',e7.2,/1x, & 'MONTH = ',i2,' HOUR = ',i3) CALL MYSTOP(1001) RETURN endif *----------------------------------------------------------------------- * solution has been found tcol = x(1) tplen = x(2) qrcol = x(3) qrwall = x(4) twall = x(5) qvwall = x(6) qvcol = x(7) qdwall = x(8) * calculate tout from energy balance tout = tplen + qvwall / ( gamma * mflow1 * aircp ) * convert kelvin to celsius tout = tout - 273.15 tcol = tcol - 273.15 tplen = tplen - 273.15 twall = twall - 273.15 troom = troom - 273.15 tamb = tamb - 273.15 tsur = tsur - 273.15 *----------------------------------------------------------------------- return end *======================================================================= * The following 14 subroutines and functions are from LAPACK, a linear * algebra package. They consist of DGESV, a general linear matrix * solver, and all of its dependents. They are (in order): * * DGESV, DGETRF, DETRS, DGETF2, DGEMM, DGER, DLASWP, DTRSM, DSCAL, * DSWAP, IDAMAX, ILAENV, LSAME, XERBLA * * DGESV is called once in the subroutine utcsolve. If the user has * another matrix solver which performs faster, they may call it from * utcsolve instead of DGESV without harming anything in the UTC system * component. If the faster solver returns the solution in [x], make * sure the lines [x]=[b] are commented out just after the call. -D.S. *======================================================================= SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * * -- LAPACK driver routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DGESV computes the solution to a real system of linear equations * A * X = B, * where A is an N-by-N matrix and X and B are N-by-NRHS matrices. * * The LU decomposition with partial pivoting and row interchanges is * used to factor A as * A = P * L * U, * where P is a permutation matrix, L is unit lower triangular, and U is * upper triangular. The factored form of A is then used to solve the * system of equations A * X = B. * * Arguments * ========= * * N (input) INTEGER * The number of linear equations, i.e., the order of the * matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the N-by-N coefficient matrix A. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (output) INTEGER array, dimension (N) * The pivot indices that define the permutation matrix P; * row i of the matrix was interchanged with row IPIV(i). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the N-by-NRHS matrix of right hand side matrix B. * On exit, if INFO = 0, the N-by-NRHS solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, so the solution could not be computed. * * ===================================================================== * * .. External Subroutines .. EXTERNAL DGETRF, DGETRS, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( N.LT.0 ) THEN INFO = -1 ELSE IF( NRHS.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -7 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGESV ', -INFO ) RETURN END IF * * Compute the LU factorization of A. * CALL DGETRF( N, N, A, LDA, IPIV, INFO ) IF( INFO.EQ.0 ) THEN * * Solve the system A*X = B, overwriting B with X. * CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, B, LDB, $ INFO ) END IF RETURN * * End of DGESV * END *======================================================================= SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DGETRF computes an LU factorization of a general M-by-N matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 3 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, IINFO, J, JB, NB * .. * .. External Subroutines .. EXTERNAL DGEMM, DGETF2, DLASWP, DTRSM, XERBLA * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRF', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * * Determine the block size for this environment. * NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN * * Use unblocked code. * CALL DGETF2( M, N, A, LDA, IPIV, INFO ) ELSE * * Use blocked code. * DO 20 J = 1, MIN( M, N ), NB JB = MIN( MIN( M, N )-J+1, NB ) * * Factor diagonal and subdiagonal blocks and test for exact * singularity. * CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO ) * * Adjust INFO and the pivot indices. * IF( INFO.EQ.0 .AND. IINFO.GT.0 ) $ INFO = IINFO + J - 1 DO 10 I = J, MIN( M, J+JB-1 ) IPIV( I ) = J - 1 + IPIV( I ) 10 CONTINUE * * Apply interchanges to columns 1:J-1. * CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) * IF( J+JB.LE.N ) THEN * * Apply interchanges to columns J+JB:N. * CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, $ IPIV, 1 ) * * Compute block row of U. * CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), $ LDA ) IF( J+JB.LE.M ) THEN * * Update trailing submatrix. * CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1, $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), $ LDA ) END IF END IF 20 CONTINUE END IF RETURN * * End of DGETRF * END *======================================================================= SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * * -- LAPACK routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. CHARACTER TRANS INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DGETRS solves a system of linear equations * A * X = B or A' * X = B * with a general N-by-N matrix A using the LU factorization computed * by DGETRF. * * Arguments * ========= * * TRANS (input) CHARACTER*1 * Specifies the form of the system of equations: * = 'N': A * X = B (No transpose) * = 'T': A'* X = B (Transpose) * = 'C': A'* X = B (Conjugate transpose = Transpose) * * N (input) INTEGER * The order of the matrix A. N >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The factors L and U from the factorization A = P*L*U * as computed by DGETRF. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (input) INTEGER array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the right hand side matrix B. * On exit, the solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL NOTRAN * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLASWP, DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 NOTRAN = LSAME( TRANS, 'N' ) IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. $ LSAME( TRANS, 'C' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( NRHS.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRS', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * IF( NOTRAN ) THEN * * Solve A * X = B. * * Apply row interchanges to the right hand sides. * CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, 1 ) * * Solve L*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, $ ONE, A, LDA, B, LDB ) * * Solve U*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, $ NRHS, ONE, A, LDA, B, LDB ) ELSE * * Solve A' * X = B. * * Solve U'*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, $ ONE, A, LDA, B, LDB ) * * Solve L'*X = B, overwriting B with X. * CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE, $ A, LDA, B, LDB ) * * Apply row interchanges to the solution vectors. * CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, -1 ) END IF * RETURN * * End of DGETRS * END *======================================================================= SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1992 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DGETF2 computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 2 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the m by n matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER J, JP * .. * .. External Functions .. INTEGER IDAMAX EXTERNAL IDAMAX * .. * .. External Subroutines .. EXTERNAL DGER, DSCAL, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETF2', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * DO 10 J = 1, MIN( M, N ) * * Find pivot and test for singularity. * JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 ) IPIV( J ) = JP IF( A( JP, J ).NE.ZERO ) THEN * * Apply the interchange to columns 1:N. * IF( JP.NE.J ) $ CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) * * Compute elements J+1:M of J-th column. * IF( J.LT.M ) $ CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) * ELSE IF( INFO.EQ.0 ) THEN * INFO = J END IF * IF( J.LT.MIN( M, N ) ) THEN * * Update trailing submatrix. * CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, $ A( J+1, J+1 ), LDA ) END IF 10 CONTINUE RETURN * * End of DGETF2 * END *======================================================================= SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * DGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of DGEMM . * END *======================================================================= SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, M, N * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of DGER . * END *======================================================================= SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. INTEGER INCX, K1, K2, LDA, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DLASWP performs a series of row interchanges on the matrix A. * One row interchange is initiated for each of rows K1 through K2 of A. * * Arguments * ========= * * N (input) INTEGER * The number of columns of the matrix A. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the matrix of column dimension N to which the row * interchanges will be applied. * On exit, the permuted matrix. * * LDA (input) INTEGER * The leading dimension of the array A. * * K1 (input) INTEGER * The first element of IPIV for which a row interchange will * be done. * * K2 (input) INTEGER * The last element of IPIV for which a row interchange will * be done. * * IPIV (input) INTEGER array, dimension (M*abs(INCX)) * The vector of pivot indices. Only the elements in positions * K1 through K2 of IPIV are accessed. * IPIV(K) = L implies rows K and L are to be interchanged. * * INCX (input) INTEGER * The increment between successive values of IPIV. If IPIV * is negative, the pivots are applied in reverse order. * * ===================================================================== * * .. Local Scalars .. INTEGER I, IP, IX * .. * .. External Subroutines .. EXTERNAL DSWAP * .. * .. Executable Statements .. * * Interchange row I with row IPIV(I) for each of rows K1 through K2. * IF( INCX.EQ.0 ) $ RETURN IF( INCX.GT.0 ) THEN IX = K1 ELSE IX = 1 + ( 1-K2 )*INCX END IF IF( INCX.EQ.1 ) THEN DO 10 I = K1, K2 IP = IPIV( I ) IF( IP.NE.I ) $ CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA ) 10 CONTINUE ELSE IF( INCX.GT.1 ) THEN DO 20 I = K1, K2 IP = IPIV( IX ) IF( IP.NE.I ) $ CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA ) IX = IX + INCX 20 CONTINUE ELSE IF( INCX.LT.0 ) THEN DO 30 I = K2, K1, -1 IP = IPIV( IX ) IF( IP.NE.I ) $ CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA ) IX = IX + INCX 30 CONTINUE END IF * RETURN * * End of DLASWP * END *======================================================================= SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of DTRSM . * END *======================================================================= subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end *======================================================================= subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end *======================================================================= integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dmax integer i,incx,ix,n c idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end *======================================================================= INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, $ N4 ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER ISPEC, N1, N2, N3, N4 * .. * * Purpose * ======= * * ILAENV is called from the LAPACK routines to choose problem-dependent * parameters for the local environment. See ISPEC for a description of * the parameters. * * This version provides a set of parameters which should give good, * but not optimal, performance on many of the currently available * computers. Users are encouraged to modify this subroutine to set * the tuning parameters for their particular machine using the option * and problem size information in the arguments. * * This routine will not function correctly if it is converted to all * lower case. Converting it to all upper case is allowed. * * Arguments * ========= * * ISPEC (input) INTEGER * Specifies the parameter to be returned as the value of * ILAENV. * = 1: the optimal blocksize; if this value is 1, an unblocked * algorithm will give the best performance. * = 2: the minimum block size for which the block routine * should be used; if the usable block size is less than * this value, an unblocked routine should be used. * = 3: the crossover point (in a block routine, for N less * than this value, an unblocked routine should be used) * = 4: the number of shifts, used in the nonsymmetric * eigenvalue routines * = 5: the minimum column dimension for blocking to be used; * rectangular blocks must have dimension at least k by m, * where k is given by ILAENV(2,...) and m by ILAENV(5,...) * = 6: the crossover point for the SVD (when reducing an m by n * matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds * this value, a QR factorization is used first to reduce * the matrix to a triangular form.) * = 7: the number of processors * = 8: the crossover point for the multishift QR and QZ methods * for nonsymmetric eigenvalue problems. * * NAME (input) CHARACTER*(*) * The name of the calling subroutine, in either upper case or * lower case. * * OPTS (input) CHARACTER*(*) * The character options to the subroutine NAME, concatenated * into a single character string. For example, UPLO = 'U', * TRANS = 'T', and DIAG = 'N' for a triangular routine would * be specified as OPTS = 'UTN'. * * N1 (input) INTEGER * N2 (input) INTEGER * N3 (input) INTEGER * N4 (input) INTEGER * Problem dimensions for the subroutine NAME; these may not all * be required. * * (ILAENV) (output) INTEGER * >= 0: the value of the parameter specified by ISPEC * < 0: if ILAENV = -k, the k-th argument had an illegal value. * * Further Details * =============== * * The following conventions have been used when calling ILAENV from the * LAPACK routines: * 1) OPTS is a concatenation of all of the character options to * subroutine NAME, in the same order that they appear in the * argument list for NAME, even if they are not used in determining * the value of the parameter specified by ISPEC. * 2) The problem dimensions N1, N2, N3, N4 are specified in the order * that they appear in the argument list for NAME. N1 is used * first, N2 second, and so on, and unused problem dimensions are * passed a value of -1. * 3) The parameter value returned by ILAENV is checked for validity in * the calling subroutine. For example, ILAENV is used to retrieve * the optimal blocksize for STRTRI as follows: * * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) * IF( NB.LE.1 ) NB = MAX( 1, N ) * * ===================================================================== * * .. Local Scalars .. LOGICAL CNAME, SNAME CHARACTER*1 C1 CHARACTER*2 C2, C4 CHARACTER*3 C3 CHARACTER*6 SUBNAM INTEGER I, IC, IZ, NB, NBMIN, NX * .. * .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL * .. * .. Executable Statements .. * GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC * * Invalid value for ISPEC * ILAENV = -1 RETURN * 100 CONTINUE * * Convert NAME to upper case if the first character is lower case. * ILAENV = 1 SUBNAM = NAME IC = ICHAR( SUBNAM( 1:1 ) ) IZ = ICHAR( 'Z' ) IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN * * ASCII character set * IF( IC.GE.97 .AND. IC.LE.122 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 10 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.97 .AND. IC.LE.122 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 10 CONTINUE END IF * ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN * * EBCDIC character set * IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN SUBNAM( 1:1 ) = CHAR( IC+64 ) DO 20 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) $ SUBNAM( I:I ) = CHAR( IC+64 ) 20 CONTINUE END IF * ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN * * Prime machines: ASCII+128 * IF( IC.GE.225 .AND. IC.LE.250 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 30 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.225 .AND. IC.LE.250 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 30 CONTINUE END IF END IF * C1 = SUBNAM( 1:1 ) SNAME = C1.EQ.'S' .OR. C1.EQ.'D' CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' IF( .NOT.( CNAME .OR. SNAME ) ) $ RETURN C2 = SUBNAM( 2:3 ) C3 = SUBNAM( 4:6 ) C4 = C3( 2:3 ) * GO TO ( 110, 200, 300 ) ISPEC * 110 CONTINUE * * ISPEC = 1: block size * * In these examples, separate code is provided for setting NB for * real and complex. We assume that NB will take the same value in * single or double precision. * NB = 1 * IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'PO' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NB = 1 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRF' ) THEN NB = 64 ELSE IF( C3.EQ.'TRD' ) THEN NB = 1 ELSE IF( C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( C2.EQ.'GB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'PB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'TR' ) THEN IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'LA' ) THEN IF( C3.EQ.'UUM' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN IF( C3.EQ.'EBZ' ) THEN NB = 1 END IF END IF ILAENV = NB RETURN * 200 CONTINUE * * ISPEC = 2: minimum block size * NBMIN = 2 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NBMIN = 8 ELSE NBMIN = 8 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF END IF ILAENV = NBMIN RETURN * 300 CONTINUE * * ISPEC = 3: crossover point * NX = 0 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( SNAME .AND. C3.EQ.'TRD' ) THEN NX = 1 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NX = 1 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF END IF ILAENV = NX RETURN * 400 CONTINUE * * ISPEC = 4: number of shifts (used by xHSEQR) * ILAENV = 6 RETURN * 500 CONTINUE * * ISPEC = 5: minimum column dimension (not used) * ILAENV = 2 RETURN * 600 CONTINUE * * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) * ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN * 700 CONTINUE * * ISPEC = 7: number of processors (not used) * ILAENV = 1 RETURN * 800 CONTINUE * * ISPEC = 8: crossover point for multishift (used by xHSEQR) * ILAENV = 50 RETURN * * End of ILAENV * END *======================================================================= LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END *======================================================================= SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * ===================================================================== * * .. Executable Statements .. * WRITE( *, FMT = 9999 )SRNAME, INFO * CALL MYSTOP(1001) RETURN * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END