1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 module Input 33 (gmat,rhside,soln_vect,wells,split,blm',mkbigvec, 34 mksparse,a_easy,a_hard,x1) where 35 36 37 import List (transpose) 38 import Matrix 39 import AbsDensematrix 40 import Utils 41 42 43 44 mkbigvec :: [[Float]] -> Vector 45 mkbigvec v = mkvector (map mkvec v) 46 47 mksparse :: [[(Int,Int,[[Float]])]] -> Matrix 48 mksparse m 49 = mkmatrix (map (map f) m) 50 where 51 f (r,c,b) = (r-1,c-1,mkblock b) 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 directory :: String 73 directory = "./Data/" 74 75 file1 :: Int -> String 76 file1 n = directory ++ "gcomp.mat" ++ (show n) ++ ".one" 77 78 file2a :: Int -> String 79 file2a n = directory ++ "gcomp.mat" ++ (show n) ++ ".twoa" 80 81 file2b :: Int -> String 82 file2b n = directory ++ "gcomp.mat" ++ (show n) ++ ".twob" 83 84 file2c :: Int -> String 85 file2c n = directory ++ "gcomp.mat" ++ (show n) ++ ".twoc" 86 87 file3 :: Int -> String 88 file3 n = directory ++ "gcomp.mat" ++ (show n) ++ ".three" 89 90 file4 :: Int -> String 91 file4 n = directory ++ "gcomp.mat" ++ (show n) ++ ".four" 92 93 file5 :: Int -> String 94 file5 n = directory ++ "gcomp.mat" ++ (show 2) ++ ".soln" 95 96 file6 :: Int -> String 97 file6 n = directory ++ "gcomp.mat" ++ (show n) ++ ".well" 98 99 file :: String 100 file = "gconem.bmtx3" -- original file 101 102 resn = 11 103 resm = 21 104 resnm = resn * resm 105 106 readmtx rtype x 107 = map (map numval) (map f y) 108 where 109 y = if (rtype == 1) then drop 1 (lines ("abc")) 110 else 111 if (rtype == 2) then init (lines ("abc")) 112 else lines ("abc") 113 114 f y = [take 15 (drop 1 y ) , 115 take 15 (drop 16 y ) , 116 take 15 (drop 31 y ) , 117 take 15 (drop 46 y ) ] 118 readsoln x 119 = map (map numval) (map f y) 120 where y = lines ("abc") 121 f y = [take 15 (drop 8 y ) , 122 take 15 (drop 24 y ) , 123 take 15 (drop 40 y ) , 124 take 15 (drop 56 y ) ] 125 126 127 myreadmtx rtype x 128 = map (map numval) (map f y) 129 where y = if (rtype == 1) then drop 1 (lines (myread x)) 130 else 131 if (rtype == 2) then init (lines (myread x)) 132 else lines (myread x) 133 f y = [take 15 (drop 1 y ) , 134 take 15 (drop 16 y ) , 135 take 15 (drop 31 y ) , 136 take 15 (drop 46 y ) ] 137 138 myread fname 139 = "abc" 140 141 142 143 readsolnvect x 144 = collect_by4 (concat z) 145 where z = map (map numval) (map f y) 146 y = lines ("abc") 147 f y = [take 21 (drop 11 y )] 148 149 collect_by4 [] = [] 150 collect_by4 (a:b:c:d:xs) 151 = [a,b,c,d]: (collect_by4 xs) 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 split m [] = [] 174 split m xs = (take m xs):(split m (drop m xs)) 175 splitmtx n m xs = split n (split m xs) 176 177 178 179 180 181 input1' n = readmtx 1 (file1 n) 182 input1 n = split 4 (input1' n) 183 184 input2a' n = readmtx 0 (file2a n) 185 input2a n = split 4 (input2a' n) 186 187 input2b' n = readmtx 0 (file2b n) 188 input2b n = split 4 (input2b' n) 189 190 input2c' n = readmtx 0 (file2c n) 191 input2c n = split 4 (input2c' n) 192 193 input3' n = readmtx 0 (file3 n) 194 input3 n = split 4 (input3' n) 195 196 input4' n = readmtx 2 (file4 n) 197 198 input5' n = readsolnvect (file5 n) 199 200 201 202 203 204 lowerband' n = take (resnm-resn) (drop resn (input1 n)) 205 lowerband n = map transpose (lowerband' n) 206 block_lowerband n = convert_lowerband (lowerband' n) 207 208 209 210 211 212 ldiagband' n = drop 1 (input2a n) 213 ldiagband n = map transpose (ldiagband n) 214 block_ldiagband n = convert_ldiagband (ldiagband' n) 215 216 217 218 219 220 mdiagband' = input2b 221 mdiagband n = map transpose (mdiagband' n) 222 block_mdiagband n = convert_mdiagband (mdiagband' n) 223 224 225 226 227 228 udiagband' n = init(input2c n) 229 udiagband n = map transpose (udiagband' n) 230 block_udiagband n = convert_udiagband (udiagband' n) 231 232 233 234 235 236 upperband' n = take (resnm-resn) (input3 n) 237 upperband n = map transpose (upperband' n) 238 block_upperband n = convert_upperband (upperband' n) 239 240 241 242 243 244 rhside n = mkbigvec (input4' n) 245 246 247 248 249 250 soln_vect n = mkbigvec (input5' n) 251 252 253 254 255 256 convert_lowerband :: [[[Float]]] -> Block_list 257 convert_lowerband b 258 = [(i+resn,i,head (drop(i-1) b) ) | i <-[1..length b]] 259 260 261 convert_ldiagband :: [[[Float]]] -> Block_list 262 convert_ldiagband b 263 = [(i+1,i,head (drop(i-1)b)) | i <-[1..length b]] 264 265 266 267 268 convert_mdiagband :: [[[Float]]] -> Block_list 269 convert_mdiagband b 270 = [(i,i,head (drop(i-1)b)) | i <-[1..length b]] 271 272 273 274 275 convert_udiagband :: [[[Float]]] -> Block_list 276 convert_udiagband b 277 = [(i,i+1,head (drop(i-1)b))| i <-[1..length b]] 278 279 280 convert_upperband :: [[[Float]]] -> Block_list 281 convert_upperband b 282 = [(i,i+resn, head (drop(i-1) b)) | i <-[1..length b]] 283 284 285 286 287 288 289 290 291 292 wells :: Int -> Int 293 wells 6 = 1 294 wells 7 = 3 295 wells 8 = 7 296 wells n = 0 297 298 gvect :: Vector 299 gvect = mkbigvec (rep resnm [1,1,1,1]) 300 301 gmat n = make_lmat (block_lowerband n) 302 (block_ldiagband n) (block_mdiagband n) (block_udiagband n) 303 (block_upperband n) 304 305 gvect0 :: Vector 306 gvect0 = mkbigvec (rep resnm [0,0,0,0]) 307 308 309 make_lmat :: Block_list -> Block_list -> Block_list -> 310 Block_list -> Block_list -> Matrix 311 make_lmat lower ldiag mdiag udiag upper 312 = mksparse [row i | i<- [1..resnm]] 313 where 314 row i = combine i lower ldiag mdiag udiag upper 315 combine i low ld md ud up = 316 if (i==1) then (take 1 md) ++ (take 1 ud) ++ (take 1 up) 317 else if (i <=resn) then diag3 ++ upper 318 else if (i >(resnm - resn)) then lower ++ diag3 319 else if (i == resnm) then 320 [last low] ++ [last ld] ++ [last md] 321 else lower ++ diag3 ++ upper 322 where 323 diag3 = (take 1 (drop (i-2) ld )) ++ 324 (take 1 (drop (i-1) md )) ++ 325 (take 1 (drop (i-1) ud )) 326 upper = take 1(drop(i-1) up) 327 lower = take 1(drop (i-resn-1) low) 328 329 330 331 332 333 334 335 x0 size = mkbigvec [ [0,0,0,0] | i<-[1..size*size]] 336 x1 size = mkbigvec [ [1,1,1,1] | i<-[1..size*size]] 337 338 a_easy size 339 = if (size > 1) then blm size 340 else mksparse [[(1,1,[[11,-1, 0, 0], 341 [-1,11,-1, 0], 342 [ 0,-1,11,-1], 343 [ 0, 0,-1,11]])]] 344 345 a_hard size 346 = if size > 1 then blm_hard size 347 else error "ummm. system size to small for this model." 348 349 350 351 off_block = [[-1,0,0,0], 352 [0,-1,0,0], 353 [0,0,-1,0], 354 [0,0,0,-1]] 355 356 d_block = [[11,-1, 0, 0], 357 [-1,11,-1, 0], 358 [ 0,-1,11,-1], 359 [ 0, 0,-1,11]] 360 361 hard_off_block = (map (map (/90)) 362 [[-1,-2,-3,-4], 363 [-4,-1,-2,-3], 364 [-3,-4,-1,-2], 365 [-2,-3,-4,-1]]) 366 367 hard_d_block = (map (map (/90)) 368 [[90,-2,-3,-4], 369 [-4,90,-2,-3], 370 [-3,-4,90,-2], 371 [-2,-3,-4,40]]) 372 373 374 375 376 377 378 379 380 381 382 383 384 385 blm :: Int -> Matrix 386 blm n = mksparse (blm' n) 387 388 blm' :: Int -> [[(Int,Int,[[Float]])]] 389 blm' n = [row i | i<- [1..(n*n)]] 390 where row i = br i n 391 392 393 br :: Int -> Int -> [(Int,Int,[[Float]])] 394 br rw n = 395 if ((rw == 1) || (rw == (n*n-n+1))) then 396 [block i| i<- [1,2,n+1,n+2]] 397 398 else if ((rw == n) || (rw == n*n)) then [block i | i<-[n-1,n,2*n-1,2*n]] 399 400 else if ((rw `mod` n) == 1) then [block i | i<-[1,2,n+1,n+2,2*n+1,2*n+2]] 401 402 else if ((rw `mod` n) == 0) then [block i | i<-[n-1,n,2*n-1,2*n,3*n-1,3*n]] 403 404 else if ((rw >1) && (rw <n)) || ((rw >(n*n-n+1)) && (rw < n*n)) then [block i | i<-[m,1+m,2+m,n+m,n+1+m,n+2+m]] 405 406 else [block i | i<-[m,1+m,2+m,n+m,n+1+m,n+2+m, 407 2*n+m,2*n+1+m,2*n+2+m]] 408 409 where 410 block i = if (rw == (i+offset)) then (rw,i+offset,d_block) 411 else (rw,i+offset,off_block) 412 m = (rw `mod` n) -1 413 offset = if (rw <= n) then 0 414 else (((rw-1) `div` n) -1) * n 415 416 417 blm_hard :: Int -> Matrix 418 blm_hard n 419 = mksparse hard 420 where 421 hard = map (map f) (blm' n) 422 f (r,c,b) = if (r/=c) then (r,c,hard_off_block) 423 else (r,c,hard_d_block) 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 include_wells m (a:as) (b:bs) (c:cs) (d:ds) 448 449 = if (as == []) then include_well m a b c d 450 else include_wells m' as bs cs ds 451 where m' = include_well m a b c d 452 453 include_well m well_pos wellg wellh welldw 454 = include_dw m'' welldw 455 where m'' = include_col_well m' well_pos wellg 456 m' = include_row_well m well_pos wellh 457 458 include_row_well m well_pos wellh 459 = m' 460 where m' = m ++ [new_row] 461 new_row = map2 f well_pos (split 4 wellh) 462 f x y = (row_num, x,[y]) 463 row_num = length m +1 464 465 466 include_col_well m well_pos wellg 467 = foldl append_row m x 468 where x = zip2 well_pos (split 4 wellg) 469 470 471 include_dw m welldw = append_row m (length m,welldw) 472 473 474 append_row m (rownum, dat) 475 = m' 476 where m' = (take (rownum-1) m) ++ 477 [appendrow] ++ 478 drop rownum m 479 appendrow = row ++ [(rownum,colnum,data')] 480 data' = split 1 dat 481 row = m !! (rownum-1) 482 colnum = length m 483 484 485 numval :: String -> Float 486 numval x = (read x)::Float 487 488 489 490 491 492 493 494 495