Return to Band.e CVS log

Up to [CSRG BSD Unix] / 41BSD / cmd / efl / efltest

Annotation of 41BSD/cmd/efl/efltest/Band.e, revision 1.1

1.1     ! root        1:   procedure bnds(n,ml,m,g,nb,b)
        !             2:  
        !             3: # to solve a*x = b, where a is a banded matrix, using gaussian
        !             4: # elimination with partial pivoting.
        !             5:  
        !             6: # mnemonic - double precision band solution of a system of
        !             7: #            linear algebraic equations.
        !             8:  
        !             9: # input -
        !            10:  
        !            11: #   n  - the order of the system.
        !            12: #   ml - the number of nonzero elements of a on and below the diagonal.
        !            13: #   m  - the total number of nonzero elements in each row of a.
        !            14: #   g  - the matrix a, with g(i,j) = a(i,i+j-ml).
        !            15: #   nb - the number of right-hand-sides b.
        !            16: #   b  - the right-hand-sides.
        !            17:  
        !            18: # output -
        !            19:  
        !            20: #   g - has been clobbered.
        !            21: #   b - the solution vectors, x.
        !            22:  
        !            23: # scratch space allocated - n*( (ml-1)*mu + 1 ) words.
        !            24:  
        !            25: # error states -
        !            26:  
        !            27: #   1 - n.lt.1.
        !            28: #   2 - ml.lt.1.
        !            29: #   3 - ml.gt.m.
        !            30: #   4 - nb.lt.1.
        !            31: #   5 - singular matrix. (recoverable)
        !            32:   
        !            33:   real g(n,m),b(n,nb)
        !            34:   integer n,ml,m,nb
        !            35:   
        !            36:   include rstack
        !            37:   
        !            38:   integer il,iint,istkgt,nerror,nerr
        !            39:   
        !            40: # check the input for errors.
        !            41:  
        !            42:   if ( n < 1 ) { seterr(" bnds - n.lt.1",14,1,2) }
        !            43:   if ( ml < 1 ) { seterr(" bnds - ml.lt.1",16,2,2) }
        !            44:   if ( ml > m ) { seterr(" bnds - ml.gt.m",15,3,2) }
        !            45:   if ( nb < 1 ) { seterr(" bnds - nb.lt.1",15,4,2) }
        !            46:  
        !            47:   enter(1)
        !            48:   
        !            49:   il = istkgt(max0(n*(ml-1),1),3) ; iint = istkgt(n,2)
        !            50:   
        !            51:   bndlu(n,ml,m,g,ws(il),is(iint))
        !            52:   
        !            53:   if ( nerror(nerr) == 0 ) { bndfb(n,ml,m,ws(il),g,is(iint),nb,b) }
        !            54:  
        !            55:   else { erroff() ; seterr(" bnds - singular matrix",23,5,1) }
        !            56:   
        !            57:   leave()
        !            58:   
        !            59:   return
        !            60:   
        !            61:   end
        !            62:   procedure bndlu(n,ml,m,g,l,int)
        !            63:  
        !            64: # to obtain the lu decomposition of a banded matrix,
        !            65: # using gaussian elimination with partial pivoting.
        !            66:  
        !            67: # mnemonic - double precision band lu decomposition.
        !            68:  
        !            69: # input -
        !            70:  
        !            71: #   n   - the order of the matrix.
        !            72: #   ml  - the number of nonzero elements of a on and below the diagonal.
        !            73: #   m   - the number of nonzero elements in each row of a.
        !            74: #   g   - the matrix a, with g(i,j) = a(i,i+j-ml).
        !            75:  
        !            76: # output -
        !            77:  
        !            78: #   l   - the lower triangular banded factor of a.
        !            79: #   g   - the upper triangular banded factor of a.
        !            80: #   int - the row pivoting used.
        !            81:  
        !            82: # scratch storage allocated - none.
        !            83:  
        !            84: # error states -
        !            85:  
        !            86: #   1 - n.lt.1.
        !            87: #   2 - ml.lt.1.
        !            88: #   3 - m.lt.ml.
        !            89: #   4 - singular matrix. (recoverable)
        !            90:   
        !            91:   real g(n,m),l(n,ml)    # l(n,ml-1).
        !            92:   integer n,ml,m,int(n)
        !            93:   
        !            94:   real x,norm,eps,r1mach
        !            95:   integer i,j,k,ll,m1,m2
        !            96:   logical sing
        !            97:   
        !            98: # check the input for errors.
        !            99:  
        !           100:   if ( n < 1 ) { seterr(" bndlu - n.lt.1",15,1,2) }
        !           101:   if ( ml < 1 ) { seterr(" bndlu - ml.lt.1",16,2,2) }
        !           102:   if ( m < ml ) { seterr(" bndlu - m.lt.ml",16,3,2) }
        !           103:  
        !           104:   entsrc(i,0)    # protect against an existing error state.
        !           105:  
        !           106:   sing = .false. ; eps = r1mach(4)
        !           107:   m1 = ml-1 ; m2 = m-ml
        !           108:   
        !           109:   ll = m1
        !           110:   for ( i = 1 , i <= min0(m1,n) , i += 1 )    # set to 0 those elements
        !           111:     {                                         # of g which are undefined.
        !           112:     do j = ml+1-i , m 
        !           113:                 { g(i,j-ll) = g(i,j) }
        !           114:     ll = ll-1
        !           115:     do j = m-ll , m 
        !           116:                 { g(i,j) = 0.0e0 }
        !           117:     }
        !           118:   
        !           119:   for ( i = 1 , i <= min0(m2,n) , i += 1 )    # zero out lower rhs wart.
        !           120:     {
        !           121:     do j = ml+i , m 
        !           122:                 { g(n+1-i,j) = 0.0e0 }
        !           123:     }
        !           124:  
        !           125:   norm = 0.0e0    # get || a || sub infinity.
        !           126:   do i = 1 , n 
        !           127:                
        !           128:     {
        !           129:     int(i) = i
        !           130:     x = 0.0e0 ; do j = 1 , m 
        !           131:                 { x += abs(g(i,j)) }
        !           132:     norm = amax1(norm,x)
        !           133:     }
        !           134:   
        !           135:   do k = 1 , n 
        !           136:                
        !           137:     {
        !           138:     x = g(k,1) ; i = k
        !           139:     
        !           140:     ll = min0(m1+k,n)
        !           141:     
        !           142:     if ( k < ll )
        !           143:       {
        !           144:       do j = k+1 , ll 
        !           145:                    # get the pivot row.
        !           146:         { if ( abs(g(j,1)) > abs(x) ) { x = g(j,1) ; i = j } }
        !           147:       }
        !           148:     
        !           149:     int(k) = i
        !           150:     
        !           151:     if ( x == 0.0e0 ) { sing = .true. ; g(k,1) = norm*eps }
        !           152:  
        !           153:     if ( ml == 1 | k == n ) { next }
        !           154:     
        !           155:     if ( i ~= k )    # need to interchange the rows.
        !           156:       {
        !           157:       do j = 1 , m 
        !           158:                 { x = g(k,j) ; g(k,j) = g(i,j) ; g(i,j) = x }
        !           159:       }
        !           160:     
        !           161:     if ( k >= ll ) { next }
        !           162:     do i = k+1 , ll 
        !           163:                
        !           164:       {
        !           165:       x = g(i,1)/g(k,1)
        !           166:       l(k,i-k) = x
        !           167:       
        !           168:       do j = 2 , m 
        !           169:                 { g(i,j-1) = g(i,j)-x*g(k,j) }
        !           170:  
        !           171:       g(i,m) = 0.0e0
        !           172:       }
        !           173:     }
        !           174:   
        !           175:   if ( sing ) { seterr(" bndlu - singular matrix",24,4,1) }
        !           176:   
        !           177:   return
        !           178:   
        !           179:   end
        !           180:   procedure bndfb(n,ml,m,l,u,int,nb,b)
        !           181:  
        !           182: # to solve l*u*x = b, where l and u result from a call to bnds.
        !           183:  
        !           184: # mnemonic - double precision band forward elimination and
        !           185: #            back-solve.
        !           186:  
        !           187: # input -
        !           188:  
        !           189: #   n   - the order of the system.
        !           190: #   ml  - the number of nonzero entries of l on and below
        !           191: #         the diagonal.
        !           192: #   m   - the number of nonzero elements of u on and above
        !           193: #         the diagonal.
        !           194: #   l   - the lower triangular banded factor.
        !           195: #   u   - the upper triangular banded factor.
        !           196: #   int - the ordering of the rows of the system, due to pivoting.
        !           197: #   nb  - the number of right-hand-sides.
        !           198: #   b   - the right-hand-sides.
        !           199:  
        !           200: # output -
        !           201:  
        !           202: #   b - the solution vectors.
        !           203:  
        !           204: # scratch space allocated - none.
        !           205:  
        !           206: # error states -
        !           207:  
        !           208: #   1 - n.lt.1.
        !           209: #   2 - ml.lt.1.
        !           210: #   3 - m.lt.ml.
        !           211: #   4 - nb.lt.1.
        !           212:  
        !           213:   real l(n,ml),u(n,m),b(n,nb)    # l(n,ml-1).
        !           214:   integer n,ml,m,int(n),nb
        !           215:  
        !           216:   integer nerror,nerr
        !           217:  
        !           218: # check the input for errors.
        !           219:  
        !           220:   if ( n < 1 ) { seterr(" bndfb - n.lt.1",15,1,2) }
        !           221:   if ( ml < 1 ) { seterr(" bndfb - ml.lt.1",16,2,2) }
        !           222:   if ( m < ml ) { seterr(" bndfb - m.lt.ml",16,3,2) }
        !           223:   if ( nb < 1 ) { seterr(" bndfb - nb.lt.1",16,4,2) }
        !           224:   
        !           225:   entsrc(nerr,0)    # protect against an existing error state.
        !           226:   
        !           227:   bndfe(n,ml,l,int,nb,b)
        !           228:                    # do the forward-elimination.
        !           229:  
        !           230:   bndbs(n,m,u,nb,b)
        !           231:                    # do the back-substitution.
        !           232:   
        !           233:   return
        !           234:   
        !           235:   end

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.