Return to Dgl.e CVS log

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

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

1.1     ! root        1:   Procedure DglsBP(Nu,Order,BC,E)
        !             2:   
        !             3: # To determine which ODE should use which boundary condition.
        !             4:   
        !             5: # Mnemonic - Double precision Galerkin's method for Linear Systems,
        !             6: #            Boundary condition Placement.
        !             7: 
        !             8: # Scratch Space Allocated -
        !             9: 
        !            10: #       S(DglsBP) <= Nu*(4*Nu+15)
        !            11: 
        !            12: # Integer words.
        !            13:   
        !            14:   Integer Nu,Order(Nu,Nu,2),BC(Nu,2,2),E(Nu,2,2)
        !            15:   
        !            16:   Integer inow,inowold,i,j,l,iMaxord,iCE,Istkgt,iPPS
        !            17:   Logical AllZero
        !            18:   
        !            19:   Struct Nodei { Integer bp,N,j,R(1) }
        !            20:   
        !            21:   Define Push  +1
        !            22:   Define Search  0
        !            23:   Define Pop  -1
        !            24: 
        !            25: # Define Node = Is(inow) -> Nodei
        !            26:   
        !            27:   Include dstack
        !            28: 
        !            29: # Check the input for errors.
        !            30: 
        !            31:   If ( Nu < 1 ) { Seterr("DglsBP - Nu.lt.1",16,1,2) }
        !            32: 
        !            33:   Do l = 1 , 2 
        !            34:     {
        !            35: 
        !            36:     Do i = 1 , Nu 
        !            37:       {
        !            38:       AllZero = .True.    # Is Order(i,.,l) = (-1,...,-1)?
        !            39:       Do j = 1 , Nu 
        !            40:         {
        !            41:         AllZero &= Order(i,j,l) == (-1)
        !            42:         If ( Order(i,j,l) < (-1) | Order(i,j,l) > 2 )
        !            43:           { Seterr("DglsBP - Order(i,j,l) not one of -1,0,1,2",41,2,2) }
        !            44:         }
        !            45:       If ( ( BC(i,1,l) ~= (-2) & BC(i,1,l) ~= 0 ) |
        !            46:            ( BC(i,2,l) ~= (-2) & BC(i,2,l) ~= 1 ) )
        !            47:         { Seterr("DglsBP - BC(i,.,l) not one of -2,0,1",36,3,2) }
        !            48:       If ( AllZero )
        !            49:         { Seterr("DglsBP - Order(i,.,l)=(-1,...,-1)",33,4,2) }
        !            50:       }
        !            51:     }
        !            52:   
        !            53:   Enter(1)
        !            54:   
        !            55:   iCE = Istkgt(Nu,2)    # Complement of E.
        !            56: 
        !            57: # Maxord(i,l) = Max over j=1,...,Nu Order(i,j,l).
        !            58: 
        !            59:   iMaxord = Istkgt(2*Nu,2)
        !            60:   Seti(2*Nu,-1,Is(iMaxord))
        !            61: 
        !            62:   Do l = 1 , 2 
        !            63:     {
        !            64:     Do i = 1 , Nu 
        !            65:       {
        !            66:       Do j = 1 , Nu 
        !            67:         {
        !            68:         Is(iMaxord+i-1+(l-1)*Nu) = Max0(Is(iMaxord+i-1+(l-1)*Nu),
        !            69:                                         Order(i,j,l))
        !            70:         }
        !            71:       }
        !            72:     }
        !            73:   
        !            74:   i = 0 ; iPPS = Push
        !            75: 
        !            76:   While ( i < 4*Nu | iPPS ~= Push )
        !            77:     {
        !            78:     
        !            79:     Switch ( iPPS )
        !            80:       {
        !            81: 
        !            82:       Case Push:    # Make a node.
        !            83: 
        !            84:         inowold = inow
        !            85:         i += 1 ; inow = Istkgt(Nu+3,2)
        !            86:         Is(inow) -> Nodei.bp = inowold
        !            87:       
        !            88: #       Get the candidates for E(i).
        !            89:       
        !            90:         D6lsBP(i,Nu,Order,BC,E,
        !            91:                Is(iMaxord),Is(iCE),Is(inow) -> Nodei.R,Is(inow) -> Nodei.N)
        !            92:       
        !            93:         Is(inow) -> Nodei.j = 0
        !            94:         iPPS = Search ; Break
        !            95: 
        !            96:       Case Search:    # Searching a node.
        !            97: 
        !            98:         Is(inow) -> Nodei.j += 1
        !            99:       
        !           100:         If ( Is(inow) -> Nodei.j > Is(inow) -> Nodei.N )    # Back-up.
        !           101:           { iPPS = Pop ; Next }
        !           102:       
        !           103:         E(i,1,1) = Is(inow) -> Nodei.R(Is(inow) -> Nodei.j)
        !           104:         iPPS = Push ; Break
        !           105: 
        !           106:       Case Pop:    # Backing up a Node.
        !           107: 
        !           108:         inow = Is(inow) -> Nodei.bp ; Istkrl(1) ; i -= 1
        !           109:         iPPS = Search ; Break
        !           110: 
        !           111:       }    # End Switch.
        !           112: 
        !           113:     If ( i == 0 )
        !           114:       { Seterr("DglsBP - Improper Boundary Conditions",37,5,1) ; Break }
        !           115: 
        !           116:     }    # End While.
        !           117: 
        !           118:   Leave()
        !           119:   
        !           120:   Return
        !           121:   
        !           122:   End
        !           123:   Procedure D6lsBP(i,Nu,Order,BC,E,
        !           124:                    Maxord,CE,R,N)
        !           125:   
        !           126:   Integer i,Nu,Order(Nu,Nu,2),BC(1),E(1),    # BC(Nu,2,2),E(Nu,2,2),
        !           127:           Maxord(Nu,2),CE(Nu),R(Nu),N    # E(i-1),R(N).
        !           128:   
        !           129:   Integer j,LR,DM,NBCs,l,ii
        !           130: 
        !           131:   If ( BC(i) < 0 ) { N = 1 ; R(N) = 0 ; Return }
        !           132: 
        !           133: # LR = 1 for left, LR = 2 for right.
        !           134: 
        !           135:   LR = 1+(i-1)/(2*Nu)
        !           136: 
        !           137: # DM = 1 for Dirichlet, DM = 2 for Mixed boundary conditions.
        !           138: 
        !           139:   DM = 1+Mod((i-1)/Nu,2)
        !           140: 
        !           141:   ii = Mod(i,Nu) ; If ( ii == 0 ) { ii = Nu }    # B(i) = B(ii,DM,LR).
        !           142:   
        !           143:   N = 0
        !           144:   
        !           145:   Do j = 1 , Nu  { CE(j) = j }    # CE = Complement of E.
        !           146:   If ( i <= 2*Nu )
        !           147:     {
        !           148:     For ( j = 1 , j < i , j += 1 )
        !           149:       {
        !           150:       If ( BC(j) >= 0 ) { CE(E(j)) = 0 }
        !           151:       }
        !           152:     }
        !           153:   Else
        !           154:     {
        !           155:     For ( j = 2*Nu+1 , j < i , j += 1 )
        !           156:       {
        !           157:       If ( BC(j) >= 0 ) { CE(E(j)) = 0 }
        !           158:       }
        !           159:     }
        !           160:   
        !           161:   Do j = 1 , Nu 
        !           162:     {
        !           163:     If ( CE(j) == 0 ) { Next }
        !           164: 
        !           165:     NBCs = 0
        !           166:     
        !           167:     For ( l = 1 , l < i , l += 1 )
        !           168:       {
        !           169:       If ( E(l) == j & BC(l) >= 0 ) { NBCs += 1 }
        !           170:       }
        !           171:     
        !           172:     If ( ( DM == 1 & Maxord(j,LR) > BC(i) ) |
        !           173:          ( DM == 2 & Order(j,ii,LR) > BC(i) ) )
        !           174:       {
        !           175:       If ( NBCs < Max0(Maxord(j,1),Maxord(j,2)) )
        !           176:         { N += 1 ; R(N) = j }
        !           177:       }
        !           178:     }
        !           179:   
        !           180:   Return
        !           181:   
        !           182:   End

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