#!/bin/sh
# This is a shell archive (produced by shar 3.49)
# To extract the files from this archive, save it to a file, remove
# everything above the "!/bin/sh" line above, and type "sh file_name".
#
# made 02/08/1991 13:35 UTC by davis@cis.ufl.edu
# Source directory /net/reef/0/davis/Distribute/SparseDyn
#
# existing files will NOT be overwritten unless -c is specified
# This format requires very little intelligence at unshar time.
# "if test", "echo", "true", and "sed" may be needed.
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 1177 -rw-r--r-- Makefile
# 4921 -rw-r--r-- README
# 364 -rw-r--r-- instructions
# 2261 -rw-r--r-- fortran-1.3.1
# 2218 -rw-r--r-- ash85
# 1398 -rw-r--r-- dwt_59
# 4807 -rw-r--r-- will199
# 1961 -rw-r--r-- will57
# 1049 -rw-r--r-- ibm32
# 4165 -rw-r--r-- cdefine.h
# 2376 -rw-r--r-- icon.h
# 705 -rw-r--r-- proc.c
# 1902 -rw-r--r-- io.c
# 3559 -rw-r--r-- dir.c
# 4491 -rw-r--r-- display.c
# 7977 -rw-r--r-- dump.c
# 3858 -rw-r--r-- mark.c
# 4100 -rw-r--r-- mattool.c
# 1582 -rw-r--r-- fdefine.h
# 6632 -rw-r--r-- global.h
# 1612 -rw-r--r-- m1.F
# 7262 -rw-r--r-- m2.F
# 2956 -rw-r--r-- m2a.F
# 2916 -rw-r--r-- m3.F
# 14875 -rw-r--r-- m4.F
# 17 -rw-r--r-- main.F
# 7313 -rw-r--r-- marko.F
# 5135 -rw-r--r-- mattool.h
# 3895 -rw-r--r-- mback.F
# 3374 -rw-r--r-- read.F
# 5008 -rw-r--r-- store.F
#
# ============= Makefile ==============
if test -f 'Makefile' -a X"$1" != X"-c"; then
echo 'x - skipping Makefile (File already exists)'
else
echo 'x - extracting Makefile (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'Makefile' &&
X#-------------------------------------------------------------------------------
X# SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X# This is freeware: you may use this tool and distribute it, as long as it's
X# free. You can contact the author at na.tdavis@na-net.ornl.gov
X#-------------------------------------------------------------------------------
XTOOLPARTS = mattool.o io.o proc.o dir.o display.o mark.o dump.o
XMARK = main.o read.o store.o m1.o m2.o m2a.o m3.o m4.o mback.o marko.o
XALL = ${TOOLPARTS} ${MARK}
XOPT = -O
XFFLAGS = ${OPT} -u
XCFLAGS = ${OPT}
XTOOLLIB = -lsuntool -lsunwindow -lpixrect
X
Xrun: d2tool
X csh instructions
X d2tool
X
Xd2tool: ${ALL} SparseMat
X f77 -o d2tool ${CFLAGS} ${ALL} ${TOOLLIB}
X# remove files no longer needed:
X rm *.o main.f read.f store.f m1.f m2.f m2a.f m3.f m4.f mback.f marko.f
X
XSparseMat:
X mkdir SparseMat
X mv will199 will57 dwt_59 ash85 ibm32 SparseMat
X compress SparseMat/*
X
X${TOOLPARTS}: mattool.h icon.h cdefine.h fdefine.h
X${MARK}: global.h fdefine.h
X
X.SUFFIXES:
X.SUFFIXES: .o .c .F
X
X.c.o:
X cc ${CFLAGS} -c $<
X
X.F.o:
X /lib/cpp $*.F | grep -v "#" | cat -s > $*.f
X f77 ${FFLAGS} -c $*.f
SHAR_EOF
true || echo 'restore of Makefile failed'
fi
# ============= README ==============
if test -f 'README' -a X"$1" != X"-c"; then
echo 'x - skipping README (File already exists)'
else
echo 'x - extracting README (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'README' &&
XSparseDynamics and the D2 algorithm (version 1.1) February 1991
X
XTim Davis
X
XComputer and Information Sciences Department,
XUniversity of Florida
X
XAbstract for NETLIB:
X--------------------------------------------------------------------------------
XSPARSDYN SparseDynamics is a sparse matrix algorithm development tool that
X runs on a Sun workstation under suntools. It allows the developer to
X observe the dynamic changes in a sparse matrix as it is operated on by
X a sparse matrix algorithm. In its current state, it only demonstrates
X an early sequential version of the D2 algorithm. SparseDynamics is a
X prototype package, so no documentation is provided on how to connect it
X to a different sparse matrix algorithm. It has been tested on a Sun
X 3/60 with 4 megabytes, a SPARCsystem 330, and a SPARCstation 1. For
X more information, refer to "A nondeterministic parallel algorithm for
X unsymmetric sparse LU factorization," T. A. Davis, and P.-C. Yew,
X SIAM Journal on Matrix Analysis and Applications (July 1990).
X
X Tim Davis, CERFACS, (European Center for Research and Advanced
X Training in Scientific Computation), Toulouse, France.
X (email: na.tdavis@na-net.ornl.gov).
X--------------------------------------------------------------------------------
X
XYou can contact the author of this program at the following address:
X
XTim Davis
XComputer and Information Sciences Department
X301 CSE, University of Florida
XGainesville, FL 32611-2024
Xphone: (904) 392-1481
Xfax: (904) 392-1220
Xemail: davis@cis.ufl.edu
X
XI can always be reached via email at na.tdavis@na-net.ornl.gov.
X
X*-------------------------------------------------------------------------------
X* SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X* This is freeware: you may use this tool and distribute it, as long as it's
X* free. You can contact the author at na.tdavis@na-net.ornl.gov
X*-------------------------------------------------------------------------------
X
XThe program consists of the following files, in a text archival form:
X
XMakefile README instructions fortran-1.3.1
Xcdefine.h fdefine.h global.h icon.h mattool.h
Xdir.c display.c dump.c io.c mark.c mattool.c proc.c
Xmain.F m1.F m2.F m2a.F m3.F m4.F marko.F mback.F read.F store.F
Xash85 ibm32 dwt_59 will199 will57
X
XThe program looks in the subdirectory "SparseMat" for all Harwell/Boeing
Xmatrices, and puts them in the "File" menu. Matrix files are all those with
X.Z tacked on the end of the filename, i.e., all compressed files in the
Xdirectory are assumed to be compressed Harwell/Boeing matrices.
XThe Makefile creates the subdirectory and places the 5 small matrices in it
X(ash85 ibm32 dwt_59 will199 will57). If you want to use more matrices, I
Xsuggest that you make a single directory somewhere with all the Harwell/Boeing
Xmatrices in it, then make a symbolic link to it from the directory containing
Xd2tool (using the command "ln -s thedirectorysomewhereelse SparseMat" in the
Xdirectory containing d2tool). The program can handle matrices of order up to
X823, a limit which is a function of the screen size (with one pixel per
Xmatrix entry).
X
XJust type "make" to compile and run, or just "d2tool" to run if you've already
Xcompiled it. If you experience difficulties compiling or running the program
Xa SPARC, take a look at the fortran-1.3.1 file and follow the instructions.
X
XThe program starts with a small, unsymmetric 6-by-6 matrix used
Xas a counter-example in the SIAM article mentioned above. Use this matrix, or
Xclick "File:" and select a different matrix. If no matrix menu appears, then
Xchange the "Dir:" to a directory with some matrices in it. Then "Run", and sit
Xback and watch, or click "Step" to single-step through the algorithm.
X
XYou can change the height of the window, but you should not change the width.
X
XButtons, toggles, etc:
X Step single-step through the algorithm, for each pivot set:
X (1) order according to number of nonzeros in row & col with
X denser rows/cols to the lower left. NOTE: this is just for
X display. The real algorithm does this symbolically.
X (2) select and hilight the pivots
X (3) permute the pivot set to the diagonal (displayed as a blank
X square).
X (4) perform the rank-k update
X
X Run run through the algorithm from the current state until the end
X
X Quit exit the program
X
X Display toggle this to display the intermediate A, L, and U matrices
X while running the algorithm
X
X Print for debugging
X
X Factor modify the pivot selection heuristic. The pivot set includes
X all pivots with MinCost <= cost <= Factor * MinCost, where
X MinCost is the pivot with lowest Markowitz cost
X If -1, then find only one pivot.
X
X Switch switch after this many pivot sets. Use
X -1 for no switch to dense. When switched,
X the dense part of LU is displayed as a big blank
X
X Dir select where the Harwell/Boeing matrices are to be found,
X defaults to SparseMat
X
X File select a matrix
SHAR_EOF
true || echo 'restore of README failed'
fi
# ============= instructions ==============
if test -f 'instructions' -a X"$1" != X"-c"; then
echo 'x - skipping instructions (File already exists)'
else
echo 'x - extracting instructions (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'instructions' &&
X#
Xclear
Xecho "***************************************************************"
Xecho "***** D2TOOL successfully compiled ***************************"
Xecho "***************************************************************"
Xecho " "
Xecho "Hit carriage return to read the instructions:"
Xecho $<
Xclear
Xmore README
Xecho "Hit carriage return to run the program:"
Xecho $<
SHAR_EOF
true || echo 'restore of instructions failed'
fi
# ============= fortran-1.3.1 ==============
if test -f 'fortran-1.3.1' -a X"$1" != X"-c"; then
echo 'x - skipping fortran-1.3.1 (File already exists)'
else
echo 'x - extracting fortran-1.3.1 (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fortran-1.3.1' &&
X
X
X Sun Fortran 1.3.1 (available only on Sparc machines)
X ====================================================
X
XNotes on path settings (set in your .cshrc):
X
X 1. Insert /usr/lang before /usr/bin and /usr/ucb in your seach PATH.
X
X 2. Insert /usr/lang/man at the start of your MANPATH.
X
X 3. If you are using OpenWindows, then insert /usr/lang/SC0.0 before
X $OPENWINHOME/lib in LD_LIBRARY_PATH.
X
X For instance, your .cshrc file might contain lines that look like this:
X
X set path=( /usr/{lang,ucb,bin} /bin /local/{bin,X11} ~/bin . )
X setenv MANPATH /usr/lang/man:/usr/man:/local/man
X setenv OPENWINHOME /usr/openwin
X setenv LD_LIBRARY_PATH /usr/lang/SC0.0:$OPENWINHOME/lib
X
X You may need to "source .cshrc" after making these changes.
X
X
X
XOther notes:
X
X * Don't mix Fortran 1.2 and fortran 1.3 libraries.
X
X * Bindings for windows are compatible with OpenWindows only through 1.0.3.
X
X
X
Xdbx bugs:
X
X * You cannot use dbx on nonreferenced dummy arguments.
X WorkAround: Be sure to reference all dummy arguments.
X
X * You cannont use dbx on structures with unions.
X
X
X
Xf77 bugs:
X
X * Fortran has problems reading from named pipes.
X
X * Redundant parentheses in I/O lists cause compiler errors.
X WorkAround: Remove the redundant parentheses.
X
X * The O and Z formats fail for characters.
X
X * Can't close /dev/null as a file.
X WorkAround: Either don't use /dev/null as a file, or don't close the unit.
X
X * The optimizer breaks function of type byte.
X WorkAround: Declare the function to be type integer instead.
X
X * I/O with records longer than 8k fails unless you compile with -Bstatic.
X
X * Complex arithmetic with -fast fails.
X
X * A common block can't be initialized in more than one program.
X WorkAround: Put both sources into the same file.
X
X * A namelist read into a type logical*1 variable always stores false.
X WorkAround: Read it into a logical*4 variable.
X
X * If you pass a constant argument to a statement functino, and if the
X definition of the statement function passes that argument in an
X expression to another function, then you sometimes get incorrect
X results.
X WorkAround: Assign the constant to a variable and pass the variable.
X
X * You can get unexpected alignment of real and integer*2 inside a union.
X WorkAround: Use integer*4.
X
X
SHAR_EOF
true || echo 'restore of fortran-1.3.1 failed'
fi
# ============= ash85 ==============
if test -f 'ash85' -a X"$1" != X"-c"; then
echo 'x - skipping ash85 (File already exists)'
else
echo 'x - extracting ash85 (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ash85' &&
X1SYMMETRIC PATTERN OF NORMAL MATRIX OF HOLLAND SURVEY. ASHKENAZI, 1974 ASH85
X 25 6 19 0 0
XPSA 85 85 304 0
X(16I5) (16I5)
X 1 6 11 14 19 22 27 32 36 42 44 49 52 55 60 66
X 69 74 78 82 85 91 94 97 102 105 110 114 117 121 126 130
X 133 137 142 147 150 153 158 163 167 170 175 179 184 188 192 195
X 199 203 207 211 215 219 223 228 232 236 239 244 246 252 257 259
X 261 263 265 268 270 272 274 276 278 280 283 285 288 290 292 294
X 296 298 300 302 304 305
X 1 2 6 7 8 2 3 8 9 10 3 4 10 4 5 10
X 11 12 5 12 23 6 7 13 14 45 7 8 14 15 16 8
X 9 16 18 9 10 11 18 19 20 10 11 11 12 20 21 22
X 12 22 23 13 44 45 14 15 45 46 47 15 16 17 31 47
X 48 16 17 18 17 18 29 30 31 18 19 28 29 19 20 27
X 28 20 21 27 21 22 24 25 26 27 22 23 24 23 24 32
X 24 25 32 33 34 25 26 34 26 27 34 36 37 27 28 37
X 39 28 29 39 29 30 39 40 30 31 40 42 43 31 43 48
X 49 32 33 64 33 34 63 64 34 35 36 64 65 35 36 38
X 65 66 36 37 38 37 38 39 38 39 66 67 69 39 40 41
X 69 70 40 41 42 71 41 70 71 42 43 50 71 72 43 49
X 50 51 44 45 52 53 54 45 46 52 57 46 47 57 58 47
X 48 58 48 49 58 59 49 51 59 60 50 51 72 73 51 60
X 61 73 52 54 57 81 53 54 55 85 54 55 56 81 55 56
X 83 84 85 56 81 82 83 57 58 80 81 58 59 80 59 60
X 61 79 80 60 61 61 62 73 76 78 79 62 73 74 76 78
X 63 64 64 65 65 66 66 67 67 68 69 68 69 69 70 70
X 71 71 72 72 73 73 74 74 75 76 75 76 76 77 78 77
X 78 78 79 79 80 80 81 81 82 82 83 83 84 84 85 85
SHAR_EOF
true || echo 'restore of ash85 failed'
fi
# ============= dwt_59 ==============
if test -f 'dwt_59' -a X"$1" != X"-c"; then
echo 'x - skipping dwt_59 (File already exists)'
else
echo 'x - extracting dwt_59 (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'dwt_59' &&
X1SYMMETRIC CONNECTION TABLE FROM DTNSRDC, WASHINGTON DWT 59
X 15 4 11 0 0
XPSA 59 59 163 0
X(16I5) (16I5) 000000000000
X 1 5 9 11 14 18 21 24 27 30 32 34 37 39 41 44
X 48 52 56 60 62 65 68 71 74 76 78 82 86 90 93 96
X 99 102 104 105 110 114 116 119 121 124 126 128 130 133 134 138
X 142 146 149 151 153 156 157 159 160 161 163 164
X 1 2 7 9 2 3 7 10 3 11 4 5 12 5 6 8
X 13 6 8 14 7 9 10 8 13 14 9 10 27 10 11 11
X 15 12 13 20 13 14 14 30 15 16 34 16 17 21 35 17
X 18 21 23 18 19 22 24 19 20 22 25 20 26 21 23 35
X 22 24 25 23 24 35 24 25 36 25 26 26 29 27 28 31
X 39 28 31 33 34 29 30 32 38 30 32 42 31 39 40 32
X 41 42 33 40 58 34 35 35 36 37 50 51 59 37 38 51
X 52 38 41 39 40 44 40 45 41 42 53 42 55 43 44 44
X 45 45 46 47 46 47 48 49 58 48 49 58 59 49 50 57
X 59 50 51 52 51 52 52 53 53 54 55 54 55 56 56 57
X 58 59 59
SHAR_EOF
true || echo 'restore of dwt_59 failed'
fi
# ============= will199 ==============
if test -f 'will199' -a X"$1" != X"-c"; then
echo 'x - skipping will199 (File already exists)'
else
echo 'x - extracting will199 (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'will199' &&
X1unsymmetric pattern of order 199 given by willoughby in reid, 1970 will199
X 57 13 44 0 0
Xpua 199 199 701 0
X(16i5) (16i5)
X 1 6 11 16 23 32 41 50 59 68 77 86 95 104 113 120
X 123 126 129 132 136 140 144 148 152 156 160 164 168 172 176 180
X 184 188 192 196 200 204 208 212 216 220 224 228 232 236 238 240
X 242 245 249 253 257 261 265 269 273 277 281 285 288 291 294 297
X 300 303 306 309 312 315 318 321 324 327 330 333 337 341 345 349
X 353 357 361 365 369 373 377 381 385 389 393 396 399 402 405 408
X 411 414 417 420 423 426 429 432 435 438 441 444 447 450 453 457
X 461 465 469 473 477 481 484 487 490 492 494 496 499 503 507 511
X 515 519 523 526 528 530 532 534 538 542 546 550 553 556 559 562
X 565 568 571 574 577 580 583 586 589 592 595 597 599 601 603 605
X 607 609 611 613 615 617 619 621 623 625 627 629 631 633 635 637
X 639 641 643 645 647 649 651 653 655 657 659 661 663 665 667 669
X 673 677 681 684 688 692 696 702
X 91 128 129 158 159 92 130 131 159 160 93 132 133 160 161 94
X 106 134 135 161 162 163 95 107 108 136 137 162 163 164 165 96
X 109 110 138 139 164 165 166 167 97 111 112 140 141 166 167 168
X 169 98 113 114 142 143 168 169 170 171 99 115 116 144 145 170
X 171 172 173 100 117 118 146 147 172 173 174 175 101 119 120 148
X 149 174 175 176 177 102 121 122 150 151 176 177 178 179 103 123
X 124 152 153 178 179 180 181 104 125 126 154 155 180 181 182 183
X 105 127 156 157 158 182 183 61 62 158 63 64 159 65 66 160
X 67 68 161 69 70 162 163 71 72 164 165 73 74 166 167 75
X 76 168 169 77 78 170 171 79 80 172 173 81 82 174 175 83
X 84 176 177 85 86 178 179 87 88 180 181 89 90 182 183 61
X 62 91 105 63 64 91 92 65 66 92 93 67 68 93 94 69
X 70 94 95 71 72 95 96 73 74 96 97 75 76 97 98 77
X 78 98 99 79 80 99 100 81 82 100 101 83 84 101 102 85
X 86 102 103 87 88 103 104 89 90 104 105 1 2 3 4 5
X 6 7 8 106 9 10 107 108 11 12 109 110 13 14 111 112
X 15 16 113 114 17 18 115 116 19 20 117 118 21 22 119 120
X 23 24 121 122 25 26 123 124 27 28 125 126 29 30 127 1
X 2 91 3 4 92 5 6 93 7 8 94 9 10 95 11 12
X 96 13 14 97 15 16 98 17 18 99 19 20 100 21 22 101
X 23 24 102 25 26 103 27 28 104 29 30 105 31 32 128 129
X 33 34 130 131 35 36 132 133 37 38 134 135 39 40 136 137
X 41 42 138 139 43 44 140 141 45 46 142 143 47 48 144 145
X 49 50 146 147 51 52 148 149 53 54 150 151 55 56 152 153
X 57 58 154 155 59 60 156 157 31 32 91 33 34 92 35 36
X 93 37 38 94 39 40 95 41 42 96 43 44 97 45 46 98
X 47 48 99 49 50 100 51 52 101 53 54 102 55 56 103 57
X 58 104 59 60 105 31 60 62 32 33 64 34 35 66 36 37
X 68 38 39 70 40 41 72 191 42 43 74 192 44 45 76 193
X 46 47 78 194 48 49 80 195 50 51 82 196 52 53 84 197
X 54 55 86 56 57 88 58 59 90 60 62 32 64 34 66 36
X 68 191 38 70 191 192 40 72 192 193 42 74 193 194 44 76
X 194 195 46 78 195 196 48 80 196 197 50 82 197 52 84 54
X 86 56 88 58 90 1 30 61 183 2 3 63 185 4 5 65
X 187 6 7 67 189 8 9 69 10 11 71 12 13 73 14 15
X 75 16 17 77 18 19 79 20 21 81 22 23 83 24 25 85
X 26 27 87 28 29 89 30 61 184 2 63 186 4 65 188 6
X 67 190 8 69 10 71 12 73 14 75 16 77 18 79 20 81
X 22 83 24 85 26 87 28 89 106 107 108 109 110 111 112 113
X 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 157
X 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
X 145 146 147 148 149 150 151 152 153 154 155 156 162 164 192 199
X 164 166 193 199 166 168 194 199 168 170 195 170 172 196 199 172
X 174 197 199 174 176 198 199 192 193 194 196 197 198
SHAR_EOF
true || echo 'restore of will199 failed'
fi
# ============= will57 ==============
if test -f 'will57' -a X"$1" != X"-c"; then
echo 'x - skipping will57 (File already exists)'
else
echo 'x - extracting will57 (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'will57' &&
X1UNSYMMETRIC PATTERN OF ORDER 57 GIVEN BY WILLOUGHBY IN REID, 1970 WILL57
X 22 4 18 0 0
XPUA 57 57 281 0
X(16I5) (16I5)
X 1 11 21 24 26 31 37 39 43 46 49 54 59 62 65 67
X 70 72 76 83 90 96 102 106 110 114 118 122 125 136 141 144
X 151 158 164 170 174 178 182 186 190 193 204 209 214 218 222 229
X 236 242 248 252 256 260 264 268 271 282
X 1 2 8 9 11 12 14 43 44 45 1 2 8 9 11 12
X 14 43 44 45 3 4 14 3 4 5 6 11 12 46 5 6
X 7 11 12 46 6 7 1 8 9 15 1 8 9 10 11 12
X 2 10 11 12 44 5 11 12 13 46 11 12 13 2 3 14
X 8 15 16 18 30 17 30 16 18 43 44 19 20 21 22 29
X 30 31 19 20 21 22 29 30 31 19 20 21 22 23 29 19
X 20 21 22 23 29 22 23 24 29 23 24 25 29 24 25 26
X 29 25 26 27 29 26 27 28 29 27 28 29 19 20 21 22
X 23 24 25 26 27 28 29 16 17 20 30 31 19 30 31 32
X 33 34 35 42 43 44 32 33 34 35 42 43 44 32 33 34
X 35 36 42 32 33 34 35 36 42 35 36 37 42 36 37 38
X 42 37 38 39 42 38 39 40 42 39 40 41 42 40 41 42
X 32 33 34 35 36 37 38 39 40 41 42 1 18 33 43 44
X 10 18 32 43 44 1 45 46 48 13 45 46 47 45 46 47
X 48 49 50 57 45 46 47 48 49 50 57 47 48 49 50 51
X 57 47 48 49 50 51 57 50 51 52 57 51 52 53 57 52
X 53 54 57 53 54 55 57 54 55 56 57 55 56 57 47 48
X 49 50 51 52 53 54 55 56 57
SHAR_EOF
true || echo 'restore of will57 failed'
fi
# ============= ibm32 ==============
if test -f 'ibm32' -a X"$1" != X"-c"; then
echo 'x - skipping ibm32 (File already exists)'
else
echo 'x - extracting ibm32 (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ibm32' &&
X1UNSYMMETRIC PATTERN ON LEAFLET ADVERTISING IBM 1971 CONFERENCE IBM32
X 11 3 8 0 0
XPUA 32 32 126 0
X(16I5) (16I5)
X 1 7 12 18 22 26 29 34 39 46 53 58 61 63 65 68
X 71 74 79 82 85 88 90 94 97 102 106 110 112 117 121 124
X 127
X 1 2 3 4 7 26 1 2 9 21 28 2 3 6 8 9
X 29 3 4 5 12 3 5 23 27 1 6 16 3 7 14 21
X 31 1 8 12 17 27 7 9 10 13 19 23 27 1 10 11
X 21 23 25 27 2 11 15 18 29 6 12 24 11 13 3 14
X 2 15 20 4 16 22 4 16 17 6 10 18 20 30 1 19
X 26 8 16 20 3 21 32 11 22 2 17 21 23 12 24 26
X 6 15 18 24 25 13 18 22 26 5 24 26 27 9 28 3
X 5 27 29 32 12 17 23 30 13 14 31 24 28 32
SHAR_EOF
true || echo 'restore of ibm32 failed'
fi
# ============= cdefine.h ==============
if test -f 'cdefine.h' -a X"$1" != X"-c"; then
echo 'x - skipping cdefine.h (File already exists)'
else
echo 'x - extracting cdefine.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cdefine.h' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#define Q(i) (i)-1
X
X#define Mmax 200000
X#define Cmax 180000
X#define Nmax 823
X#define Nil 0
X#define Bl 32
X#define Bl1 (Bl+1)
X#define Rowbl (Bl+1)
X#define Colbl (Bl+2)
X
X#define Im(i) (global_.im [Q(i)])
X#define Xm(i) (global_.xm [Q(i)])
X#define MemHead (global_.memhead)
X#define MemTail (global_.memtail)
X#define MemLock (global_.memlock)
X
X#define Cm(i) (global_.cm [Q(i)])
X#define MemCHead (global_.memchead)
X#define MemCTail (global_.memctail)
X#define MemCLock (global_.memclock)
X
X#define Size (global_.size)
X#define Nonzeros (global_.nonzeros)
X#define Name (global_.name)
X
X#define RowLen(i) (global_.rowlen [Q(i)])
X#define RowHead(i) (global_.rowhead [Q(i)])
X
X#define ColLen(i) (global_.collen [Q(i)])
X#define ColHead(i) (global_.colhead [Q(i)])
X
X#define Lp(i) (global_.lp [Q(i)])
X#define Up(i) (global_.up [Q(i)])
X#define Llen(i) (global_.llen [Q(i)])
X#define Ulen(i) (global_.ulen [Q(i)])
X
X#define Pivrow(i) (global_.pivrow [Q(i)])
X#define Ipivrow(i) (global_.ipivrow [Q(i)])
X#define Pivcol(i) (global_.pivcol [Q(i)])
X#define Ipivcol(i) (global_.ipivcol [Q(i)])
X
X#define Iw1(i) (global_.iw1 [Q(i)])
X#define Iw2(i) (global_.iw2 [Q(i)])
X#define Iw3(i) (global_.iw3 [Q(i)])
X#define Iw4(i) (global_.iw4 [Q(i)])
X#define Wi(i) (global_.wi [Q(i)])
X#define W(i) (global_.w [Q(i)])
X#define Rhs(i) (global_.rhs [Q(i)])
X#define Xt(i) (global_.xt [Q(i)])
X#define Error (global_.error)
X#define Lnz (global_.lnz)
X#define Unz (global_.unz)
X#define Prows(i) (global_.prows [Q(i)])
X#define Pcols(i) (global_.pcols [Q(i)])
X#define Stagesize(i) (global_.stagesize [Q(i)])
X#define Stage (global_.stage)
X#define NStages (global_.nstages)
X#define Printflag (global_.printflag)
X#define Debug (global_.debug)
X#define Sdense (global_.sdense)
X#define Idense (global_.idense)
X#define Densep (global_.densep)
X#define Densesize (global_.densesize)
X#define Sparsesize (global_.sparsesize)
X#define Stage (global_.stage)
X#define NPivots (global_.npivots)
X#define Pfirst (global_.pfirst)
X#define Plast (global_.plast)
X#define Redrow(i) (global_.redrow [Q(i)])
X#define Redcol(i) (global_.redcol [Q(i)])
X#define NRedrow (global_.nredrow)
X#define NRedcol (global_.nredcol)
X#define Lpointer(i) (global_.pcols [Q(i)])
X
X#define MRowhead(i) (global_.mrowhead [Q(i)])
X#define MRownext(i) (global_.mrownext [Q(i)])
X#define MRowlast(i) (global_.mrowlast [Q(i)])
X#define MColhead(i) (global_.mcolhead [Q(i)])
X#define MColnext(i) (global_.mcolnext [Q(i)])
X#define MCollast(i) (global_.mcollast [Q(i)])
X#define MRowmin (global_.mrowmin)
X#define MColmin (global_.mcolmin)
X#define MRowflag(i) (global_.mrowflag [Q(i)])
X#define MColflag(i) (global_.mcolflag [Q(i)])
X#define MRowmax(i) (global_.mrowmax [Q(i)])
X#define MFactor (global_.mfactor)
X
X/* #define Z(i,j) (global_.z [Q(j)][Q(i)]) */
X
Xtypedef struct global_data
X{
X double
X/* z [Nmax][Nmax], */
X xm [Mmax],
X rhs [Nmax], xt [Nmax], w [Nmax+2], mrowmax [Nmax], mfactor, error,
X sdense
X ;
X int
X im [Mmax], memhead, memtail, memlock,
X cm [Cmax], memchead, memctail, memclock,
X size, nonzeros,
X rowlen [Nmax], rowhead [Nmax],
X collen [Nmax], colhead [Nmax],
X lp [Nmax], up [Nmax], llen [Nmax], ulen [Nmax],
X pivrow [Nmax], ipivrow [Nmax],
X pivcol [Nmax], ipivcol [Nmax],
X mrowhead [Nmax], mrownext [Nmax], mrowlast [Nmax],
X mcolhead [Nmax], mcolnext [Nmax], mcollast [Nmax],
X mrowmin, mcolmin, mrowflag [Nmax], mcolflag [Nmax],
X prows [Nmax], pcols [Nmax], stagesize [Nmax],
X iw1 [Nmax+2], iw2 [Nmax+2], iw3 [Nmax+2], iw4 [Nmax+2], wi [Nmax+2],
X stage, nstages, printflag, npivots, pfirst, plast, debug,
X idense, densesize, densep, sparsesize,
X redrow [Nmax], redcol [Nmax], nredrow, nredcol,
X lnz, unz
X ;
X char
X name [10]
X ;
X}
Xglobal_data ; extern global_data global_ ;
X
X#include "fdefine.h"
X#define idint(x) ((int) x)
SHAR_EOF
true || echo 'restore of cdefine.h failed'
fi
# ============= icon.h ==============
if test -f 'icon.h' -a X"$1" != X"-c"; then
echo 'x - skipping icon.h (File already exists)'
else
echo 'x - extracting icon.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'icon.h' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov.
X This file was modified from an icon by Allan Tuchman.
X------------------------------------------------------------------------------*/
X/* Format_version=1, Width=64, Height=64, Depth=1, Valid_bits_per_item=16
X */
X 0x0000,0x0000,0x0000,0x0000,0x7FF0,0x0000,0x0000,0x0FFE,
X 0x7FF0,0x0000,0x0000,0x0FFE,0x6000,0x0000,0x0000,0x0006,
X 0x6000,0x0000,0x0000,0x0006,0x6555,0x554F,0xFFFF,0xFFE6,
X 0x62AA,0xAAA8,0x0000,0x0026,0x6155,0x5548,0x0000,0x0026,
X 0x60AA,0xAAA8,0x0000,0x0026,0x6455,0x5548,0x0000,0x0026,
X 0x662A,0xAAA8,0x0000,0x0026,0x6515,0x5548,0x0000,0x0026,
X 0x648A,0xAAA8,0x0000,0x0026,0x6445,0x5548,0x0000,0x0026,
X 0x6422,0xAAA8,0x0000,0x0026,0x6411,0x5548,0x0000,0x0026,
X 0x6408,0xAAA8,0x0000,0x0026,0x6404,0x5548,0x0000,0x0026,
X 0x6402,0x2AA8,0x0000,0x0026,0x6401,0x1548,0x0000,0x0026,
X 0x6400,0x8AA8,0x0000,0x0026,0x6400,0x4548,0x0000,0x0026,
X 0x6400,0x22A8,0x0000,0x0026,0x6400,0x1148,0x0000,0x0026,
X 0x6400,0x08A8,0x0000,0x0026,0x6400,0x0448,0x0000,0x0026,
X 0x6400,0x0228,0x0000,0x0026,0x6400,0x0108,0x0000,0x0026,
X 0x6400,0x0088,0x0000,0x0026,0x6400,0x0048,0x0000,0x0026,
X 0x6400,0x0028,0x0000,0x0026,0x6400,0x0028,0x0000,0x0026,
X 0x6403,0x8028,0x0000,0x0026,0x6407,0xC028,0x0000,0x0026,
X 0x6406,0xC028,0x0000,0x0026,0x640E,0xE028,0x0000,0x0026,
X 0x640E,0xE028,0x0000,0x0026,0x640E,0xE028,0x0000,0x0026,
X 0x640E,0xE028,0x0000,0x0026,0x640E,0xE028,0x0000,0x0026,
X 0x6406,0xC028,0x0000,0x0026,0x6407,0xC028,0x0000,0x0026,
X 0x6403,0x8028,0x0000,0x0026,0x6400,0x0028,0x0000,0x0026,
X 0x6400,0x0028,0x0000,0x0026,0x6400,0x0028,0x0000,0x0026,
X 0x6400,0x0028,0x0000,0x0026,0x67FF,0xFFFF,0xFFFF,0xFFE6,
X 0x6000,0x0000,0x0000,0x0006,0x6000,0x0000,0x0000,0x0006,
X 0x7FF0,0x0000,0x0000,0x0FFE,0x7FF0,0x0000,0x0000,0x0FFE,
X 0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,
X 0x0000,0x0181,0x8000,0x03C0,0x0000,0x0181,0x8000,0x00C0,
X 0x0EC3,0xC7E7,0xE3C3,0xC0C0,0x0FE6,0x6181,0x8666,0x60C0,
X 0x0D60,0x6181,0x8666,0x60C0,0x0D63,0xE181,0x8666,0x60C0,
X 0x0D66,0x6181,0x8666,0x60C0,0x0D66,0x61A1,0xA666,0x60C0,
X 0x0D63,0xE0C0,0xC3C3,0xC0C0,0x0000,0x0000,0x0000,0x0000
SHAR_EOF
true || echo 'restore of icon.h failed'
fi
# ============= proc.c ==============
if test -f 'proc.c' -a X"$1" != X"-c"; then
echo 'x - skipping proc.c (File already exists)'
else
echo 'x - extracting proc.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'proc.c' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include "mattool.h"
X
Xprint_proc (item, value, event)
XPanel_item item ;
Xunsigned int value ;
XEvent *event ;
X{
X Printflag = value ? TRUE : FALSE ;
X}
X
Xfast_proc (item, value, event)
XPanel_item item ;
Xunsigned int value ;
XEvent *event ;
X{
X display = value ? TRUE : FALSE ;
X}
X
Xquit_proc ()
X{
X window_destroy (frame);
X}
SHAR_EOF
true || echo 'restore of proc.c failed'
fi
# ============= io.c ==============
if test -f 'io.c' -a X"$1" != X"-c"; then
echo 'x - skipping io.c (File already exists)'
else
echo 'x - extracting io.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'io.c' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include "mattool.h"
X
Xload_proc ()
X{
X if (init_file () && read_harwell ())
X {
X system ("rm matrix.in") ;
X MFactor = dget_value (factor_item, 2.0) ;
X mstep = 1 ;
X printf ("PMark:%8.2f %8s %s\n", MFactor, filename, path) ;
X status = Loaded ;
X }
X else
X {
X printf ("Invalid matrix: %s\n", path) ;
X status = None ;
X }
X
X}
X
Xinit_file ()
X{
X char command [MAXPATHLEN] ;
X
X if (*pathname == '\0')
X {
X printf ("invalid directory specification\n") ;
X window_set (control_panel, PANEL_CARET_ITEM, dir_item, 0);
X return (FALSE);
X }
X if (*filename == '\0')
X {
X printf ("invalid filename specification\n") ;
X return (FALSE);
X }
X sprintf (path, "%s/%s", pathname, filename);
X sprintf (command, "zcat %s > matrix.in\n", path) ;
X if (system (command) != 0)
X {
X printf ("Can't open %s\n", filename) ;
X return (FALSE) ;
X }
X return (TRUE) ;
X}
X
Xread_harwell ()
X{
X int j, k, i, p, row, col ;
X
X readmatrix_ () ;
X if (Size <= 0) return (0) ;
X scale_display () ;
X display_original () ;
X storematrix_ () ;
X return (1) ;
X}
X
Xdouble dget_value (item, dfault)
XPanel_item item ;
Xdouble dfault ;
X{
X char *p, q [20] ;
X double value ;
X
X p = (char *) panel_get_value (item) ;
X if (p)
X {
X if (*p == '*') value = -1 ;
X else sscanf (p, "%lg", &value) ;
X }
X else value = dfault ;
X
X if (value < 0.0) panel_set_value (item, "*") ;
X else panel_set_value (item, sprintf (q, "%g", value)) ;
X
X return (value) ;
X}
X
SHAR_EOF
true || echo 'restore of io.c failed'
fi
# ============= dir.c ==============
if test -f 'dir.c' -a X"$1" != X"-c"; then
echo 'x - skipping dir.c (File already exists)'
else
echo 'x - extracting dir.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'dir.c' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include "mattool.h"
X
Xdir_proc (item, event)
XPanel_item item ;
XEvent *event ;
X{
X char *dir ;
X
X dir = (char *) panel_get_value (item) ;
X if (*dir == NULL || NCompare (dir, " ", 1))
X {
X *pathname = '\0' ;
X printf ("No directory specified\n") ;
X return ;
X }
X if (!Compare (pathname, dir))
X {
X strcpy (pathname, dir) ;
X dir_setup () ;
X }
X}
X
Xdir_setup ()
X{
X struct direct *dp ;
X int sparse, matsize, m, i, j, k, slot ;
X char command [MAXPATHLEN], s [100], *z, *dot ;
X FILE *matrix_file ;
X Menu_item item ;
X
X printf ("directory: %s\n", pathname) ;
X
X if ((cwd = opendir (pathname)) == NULL) /* open directory */
X {
X printf ("Can't open directory\n") ;
X return ;
X }
X
X if (file_menu) menu_destroy (file_menu) ;
X file_menu = menu_create (MENU_NCOLS, 6, 0) ;
X entries = 0 ;
X
X while ((dp = readdir (cwd)) != NULL) /* read next directory entry */
X {
X z = rindex (dp->d_name, 'Z') ; /* find .Z in file name */
X dot = rindex (dp->d_name, '.') ;
X
X /* sparse if ".Z" are last two characters of file name */
X sparse = (z != NULL) && (dot != NULL) && (z == dot + 1) ;
X
X if (sparse && entries < Maxfiles-1)
X {
X *dot = '\0' ; /* eliminate trailing ".Z" */
X
X /* create pipe from "zcat filename" command */
X sprintf (command, "zcat %s/%s\n", pathname, dp->d_name) ;
X matsize = 0 ;
X if ((matrix_file = popen (command, "r")) == NULL)
X printf ("Can't open %s\n", dp->d_name) ;
X else
X {
X /* read size of matrix */
X fgets (s, Line_len, matrix_file) ;
X fgets (s, Line_len, matrix_file) ;
X fgets (s, Line_len, matrix_file) ;
X sscanf (s, "%3s %14d", type, &matsize) ;
X
X if (matsize <= 0)
X printf ("\n Can't read %s\n", dp->d_name) ;
X else if (matsize <= Nmax)
X {
X entries++ ;
X strcpy (files [entries], dp->d_name) ;
X msize [entries] = matsize ;
X }
X }
X
X /* close the file */
X pclose (matrix_file) ;
X }
X }
X closedir (cwd) ; /* close the directory */
X
X slot = 0 ; /* sort by matrix size */
X for (i = 1 ; i <= entries ; i++)
X {
X k = i ;
X for (j = i+1 ; j <= entries ; j++)
X if (msize [j] < msize [k] ||
X (msize [j] == msize [k] && strcmp (files [j], files [k]) < 0))
X k = j ;
X
X if (i != k)
X {
X m = msize [i] ;
X msize [i] = msize [k] ;
X msize [k] = m ;
X strcpy (s, files [i]) ;
X strcpy (files [i], files [k]) ;
X strcpy (files [k], s) ;
X }
X printf ("%s %d\n", files [i], msize [i]) ;
X
X menu_set (file_menu, MENU_STRING_ITEM, files [i], ++slot, 0) ;
X sprintf (smsize [i], "%d", msize [i]) ;
X menu_set (file_menu, MENU_STRING_ITEM, smsize [i], ++slot, 0) ;
X item = menu_get (file_menu, MENU_NTH_ITEM, slot) ;
X menu_set (item, MENU_INACTIVE, TRUE, 0) ;
X }
X Clear_msg (file_item) ;
X}
X
Xfile_proc (item, event)
XPanel_item item ;
XEvent *event ;
X{
X int i = ((int) menu_show (file_menu, control_panel, event, 0) + 1) / 2 ;
X if (i < 1) return ;
X strcpy (filename, files [i]) ;
X Set_msg (file_item, filename) ;
X load_proc () ;
X}
SHAR_EOF
true || echo 'restore of dir.c failed'
fi
# ============= display.c ==============
if test -f 'display.c' -a X"$1" != X"-c"; then
echo 'x - skipping display.c (File already exists)'
else
echo 'x - extracting display.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'display.c' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include "mattool.h"
X
Xdisplay_either ()
X{
X if (status == Loaded) display_original () ;
X else if (status == Done || status == Running) display_ALU () ;
X}
X
Xscale_display ()
X{
X pw = canvas_pixwin (canvas);
X wc = (int) window_get (canvas, CANVAS_WIDTH) ;
X wr = (int) window_get (canvas, CANVAS_HEIGHT) ;
X printf ("Window size: wc=%d wr=%d\n", wc, wr) ;
X
X n = Size ;
X wsize = Min (wc, wr) - 4 ;
X ratio = wsize / n ;
X wsize = (ratio * n) + 4 ;
X elsize = (ratio <= 2) ? ratio : ratio - 1 ;
X /*
X printf ("Wsize = %d\n", wsize);
X printf ("elsize: %d ratio: %d\n", elsize, ratio) ;
X */
X if (ratio == 0) { printf ("Matrix too large\n") ; return (FALSE) ; }
X
X clear_canvas () ;
X
X pw_vector (pw, 0, wsize-1, wsize-1, wsize-1, PIX_SRC, 1) ;
X pw_vector (pw, wsize-1, wsize-1, wsize-1, 0, PIX_SRC, 1) ;
X return (TRUE) ;
X}
X
Xclear_canvas ()
X{
X Rect *rect ;
X
X rect = (Rect *) window_get (canvas, WIN_RECT) ;
X pw_writebackground (pw, 0, 0, rect->r_width, rect->r_height, PIX_SRC) ;
X}
X
Xdisplay_original ()
X{
X int i, row, col, k ;
X
X pw_batch_on (pw) ;
X display_clear () ;
X for (k = 1 ; k <= Nonzeros ; k++) Set_element (A_row (k), A_col (k)) ;
X pw_batch_off (pw) ;
X}
X
Xdisplay_ALU ()
X{
X pw_batch_on (pw) ;
X display_clear () ;
X display_A_rows () ;
X if (NStages > 0)
X {
X display_L () ;
X display_U () ;
X display_box () ;
X if (Densep == 0 && mstep != 1) show_updates () ;
X }
X pw_batch_off (pw) ;
X}
X
Xdisplay_box ()
X{
X int i, h, diag ;
X
X diag = 1 ;
X for (i = 1 ; i <= NStages ; i++)
X {
X h = ratio * Stagesize (i) ;
X Clear (Scrx (diag), Scry (diag), h, h) ;
X Set (Scrx (diag), Scry (diag), h, 1) ;
X Set (Scrx (diag), Scry (diag), 1, h) ;
X Set (Scrx (diag), Scry (diag)+h-1, h, 1) ;
X Set (Scrx (diag)+h-1, Scry (diag), 1, h) ;
X diag += Stagesize (i) ;
X }
X}
X
Xdisplay_A ()
X{
X int row, col, j, k, p, len, i ;
X
X if (Densep != 0) return ;
X
X for (i = Stage+1 ; i <= Sparsesize ; i++)
X {
X k = 1 ;
X col = Pivcol (i) ;
X p = ColHead (col) ;
X len = ColLen (col) ;
X for (j = 1 ; j <= len ; j++)
X {
X row = Col_index (p, k) ;
X Set_element (Ipivrow (row), Ipivcol (col)) ;
X k++ ;
X if (k > Bl) { k = 1 ; p = Col_next (p) ; }
X }
X }
X}
X
Xdisplay_A_rows ()
X{
X int row, col, j, k, p, len, i ;
X
X if (Densep != 0) return ;
X
X for (i = Stage+1 ; i <= Sparsesize ; i++)
X {
X k = 1 ;
X row = Pivrow (i) ;
X p = RowHead (row) ;
X len = RowLen (row) ;
X for (j = 1 ; j <= len ; j++)
X {
X col = Row_index (p, k) ;
X Set_element (Ipivrow (row), Ipivcol (col)) ;
X k++ ;
X if (k > Bl) { k = 1 ; p = Row_next (p) ; }
X }
X }
X}
X
Xdisplay_L ()
X{
X int row, col, p, len, j, s ;
X
X s = (Densep == 0) ? Stage : Sparsesize ;
X
X for (col = 1 ; col <= s ; col++)
X {
X p = Lp (col) ;
X len = Llen (col) ;
X for (j = 1 ; j <= len ; j++)
X {
X row = L_index (p, j) ;
X Set_element (Ipivrow (row), col) ;
X }
X }
X}
X
Xdisplay_U ()
X{
X int row, col, p, len, j, s ;
X
X s = (Densep == 0) ? Stage : Sparsesize ;
X
X for (row = 1 ; row <= s ; row++)
X {
X p = Up (row) ;
X len = Ulen (row) ;
X for (j = 1 ; j <= len ; j++)
X {
X col = U_index (p, j) ;
X Set_element (row, Ipivcol (col)) ;
X }
X }
X}
X
Xdisplay_pivots ()
X{
X int i, row, col, hisize = (ratio > 2) ? ratio : 3 ;
X
X pw_batch_on (pw) ;
X for (i = 1 ; i <= NPivots ; i++)
X {
X col = Ipivcol (Pcols (i)) ;
X row = Ipivrow (Prows (i)) ;
X Set (Scrx (col), Scry (row), hisize, hisize) ;
X Clear (Scrx (col)+1, Scry (row)+1, hisize-2, hisize-2) ;
X }
X pw_batch_off (pw) ;
X}
X
Xshow_updates ()
X{
X int rows, cols ;
X
X cols = Plast + NRedcol ;
X rows = Plast + NRedrow ;
X Set (Scrx (cols + 1) - 1, Scry (1), 1, ratio * rows) ;
X Set (Scrx (1), Scry (rows + 1) - 1, ratio * cols, 1) ;
X Set (Scrx (Plast+1)-1, Scry (Plast+1)-1, ratio * (Size - Plast), 1) ;
X Set (Scrx (Plast+1)-1, Scry (Plast+1)-1, 1, ratio * (Size - Plast)) ;
X}
X
Xdisplay_clear ()
X{
X Clear (Scrx (1), Scry (1), (Size*ratio)+1, (Size*ratio)+1) ;
X}
SHAR_EOF
true || echo 'restore of display.c failed'
fi
# ============= dump.c ==============
if test -f 'dump.c' -a X"$1" != X"-c"; then
echo 'x - skipping dump.c (File already exists)'
else
echo 'x - extracting dump.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'dump.c' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include "mattool.h"
X
Xdumpmatrix_ ()
X{
X int j, k, p, row, col, len, prow, pcol ;
X
X dumpu_ () ;
X dumpl_ () ;
X
X printf ("\n\nSparse Columns of A:") ;
X for (col = Stage+1 ; col <= Sparsesize ; col++)
X {
X k = 1 ;
X pcol = Pivcol (col) ;
X p = ColHead (pcol) ;
X len = ColLen (pcol) ;
X
X printf ("\n%3d: pcol: %d len: %d p: %d", col, pcol, len, p) ;
X if (len > Bl) printf (" LARGE *****") ;
X printf ("\n ") ;
X
X for (j = 1 ; j <= len ; j++)
X {
X prow = Ipivrow (Col_index (p, k)) ;
X printf (" %3d", prow) ;
X k++ ;
X if (k > Bl) { k = 1 ; p = Col_next (p) ; }
X }
X }
X
X printf ("\n\nSparse Rows of A: stage %d", Stage) ;
X for (row = Stage+1 ; row <= Sparsesize ; row++)
X {
X k = 1 ;
X prow = Pivrow (row) ;
X p = RowHead (prow) ;
X len = RowLen (prow) ;
X
X printf ("\n%3d: prow: %d len: %d p: %d", row, prow, len, p) ;
X if (len > Bl) printf (" LARGE *****") ;
X printf ("\n ") ;
X
X for (j = 1 ; j <= len ; j++)
X {
X pcol = Ipivcol (Row_index (p, k)) ;
X printf ("(%3d %3d %g)", pcol, Row_index (p, k), Row_value (p, k)) ;
X k++ ;
X if (k > Bl) { k = 1 ; p = Row_next (p) ; }
X }
X }
X printf ("\n") ;
X}
X
Xdumpnzlist_ ()
X{
X int row, col, i, j, nz, first ;
X
X printf ("\nRow nz: %d, size %d", MRowmin, Size) ;
X first = TRUE ;
X for (nz = 1 ; nz <= Size ; nz++)
X {
X row = MRowhead (nz) ;
X if (row != 0)
X {
X if (first && (nz != MRowmin))
X {
X printf ("ERROR: MRowmin\n") ;
X crwait_ () ;
X }
X first = FALSE ;
X printf ("\n%3d: ", nz) ;
X while (row != 0)
X {
X printf (" (%3d %3d %3d)", MRowlast (row), row, MRownext (row)) ;
X/* printf (" %3d", row) ; */
X if (MRowflag (row) == 1)
X {
X printf ("\nERROR: %3d flagged\n", row) ;
X crwait_ () ;
X }
X if (RowLen (row) != nz)
X {
X printf ("\nERROR: RowLen(%d)=%d\n", row, RowLen (row)) ;
X crwait_ () ;
X }
X row = MRownext (row) ;
X }
X }
X }
X
X printf ("\nCol nz: %d", MColmin) ;
X for (nz = 1 ; nz <= Size ; nz++)
X {
X col = MColhead (nz) ;
X if (col != 0)
X {
X if (first && (nz != MColmin))
X {
X printf ("ERROR: MColmin\n") ;
X crwait_ () ;
X }
X first = FALSE ;
X printf ("\n%3d: ", nz) ;
X while (col != 0)
X {
X printf (" (%3d %3d %3d)", MCollast (col), col, MColnext (col)) ;
X/* printf (" %3d", col) ; */
X if (MColflag (col) == 1)
X {
X printf ("\nERROR: %3d flagged\n", col) ;
X crwait_ () ;
X }
X if (ColLen (col) != nz)
X {
X printf ("\nERROR: ColLen(%d)=%d\n", col, ColLen (col)) ;
X crwait_ () ;
X }
X col = MColnext (col) ;
X }
X }
X }
X printf ("\n") ;
X}
X
Xdumpu_ ()
X{
X int row, p, len, j, s ;
X
X s = (Densep == 0) ? Stage : Sparsesize ;
X printf ("\n\nSparse Rows of U: s=%d", s) ;
X for (row = 1 ; row <= s ; row++)
X {
X p = Up (row) ;
X len = Ulen (row) ;
X printf ("\n%3d: p: %7d len: %3d\n ", row, p, len) ;
X for (j = 1 ; j <= len ; j++)
X printf ("(%3d %3d) ", U_index (p, j), Ipivcol (U_index (p, j))) ;
X }
X}
X
Xdumpl_ ()
X{
X int col, p, len, j, s ;
X
X s = (Densep == 0) ? Stage : Sparsesize ;
X printf ("\n\nSparse Columns of L: s=%d", s) ;
X for (col = 1 ; col <= s ; col++)
X {
X p = Lp (col) ;
X len = Llen (col) ;
X printf ("\n%3d: p: %7d len: %3d\n ", col, p, len) ;
X for (j = 1 ; j <= len ; j++) printf (" %3d", L_index (p, j)) ;
X }
X}
X
Xdumpllist_ ()
X{
X int row, p, len, j, i ;
X
X printf ("\n\nLinked rows in the current Ls:\n") ;
X for (i = 1 ; i <= NRedrow ; i++)
X {
X row = Redrow (i) ;
X printf ("%3d %3d:\n", row, Ipivrow (row)) ;
X for (p = Lpointer (row) ; p != Nil ; p = L_next (p))
X printf (" ( %7d %4d)", p, (int) L_col (p)) ;
X printf ("\n") ;
X }
X}
X
Xdumprow_ (row)
Xint *row ;
X{
X int k, p, len, col, j ;
X double amax = 0.0 ;
X
X k = 1 ;
X p = RowHead (*row) ;
X len = RowLen (*row) ;
X
X printf ("\nrow: (%4d %4d) len: %4d p: %8d\n", *row, Ipivrow (*row), len, p) ;
X
X for (j = 1 ; j <= len ; j++)
X {
X col = Row_index (p, k) ;
X printf (" %4d %4d %8d (%4d %4d) %14.6g\n",
X j, k, p, col, Ipivcol (col), Row_value (p, k)) ;
X amax = Max (amax, Abs (Row_value (p, k))) ;
X k++ ;
X if (k > Bl) { k = 1 ; p = Row_next (p) ; }
X }
X printf (" maximum value: %14.6g\n", amax) ;
X}
X
Xdumpcol_ (col)
Xint *col ;
X{
X int k, p, len, row, j ;
X
X k = 1 ;
X p = ColHead (*col) ;
X len = ColLen (*col) ;
X
X printf ("\ncol: (%4d %4d) len: %4d p: %8d\n", *col, Ipivcol (*col), len, p) ;
X
X for (j = 1 ; j <= len ; j++)
X {
X row = Col_index (p, k) ;
X printf (" %4d %4d %8d (%4d %4d)\n", j, k, p, row, Ipivrow (row)) ;
X k++ ;
X if (k > Bl) { k = 1 ; p = Col_next (p) ; }
X }
X}
X
Xcheckrowcol_ ()
X{
X int k, p, len, row, col, j, i ;
X
X printf ("start checkrowcol:\n") ;
X for (i = 1 ; i <= Size ; i++) Iw1 (i) = 0 ;
X
X for (i = Stage+1 ; i <= Size ; i++)
X {
X k = 1 ;
X row = Pivrow (i) ;
X p = RowHead (row) ;
X len = RowLen (row) ;
X
X for (j = 1 ; j <= len ; j++)
X {
X Iw1 (Row_index (p, k))++ ;
X k++ ;
X if (k > Bl) { k = 1 ; p = Row_next (p) ; }
X }
X }
X
X for (i = Stage+1 ; i <= Size ; i++)
X {
X col = Pivcol (i) ;
X if (Iw1 (col) != ColLen (col))
X printf ("\nERROR: row / col %d mismatch: by row %d, by col %d\n",
X col, Iw1 (col), ColLen (col)) ;
X }
X printf ("done checkrowcol\n") ;
X}
X
Xdumpperm_ ()
X{
X int i ;
X
X printf ("permutation matrices:\n") ;
X for (i = 1 ; i <= Size ; i++)
X printf ("%3d: %3d %3d %3d %3d\n", i, Pivrow (i), Ipivrow (i),
X Pivcol (i), Ipivcol (i)) ;
X}
X
Xdumprhs_ ()
X{
X int i ;
X
X printf ("Right hand side (unpermuted / row-permuted)\n") ;
X for (i = 1 ; i <= Size ; i++)
X printf ("%3d: %20.10g %3d %20.10g\n",
X i, Rhs (i), Pivrow (i), Rhs (Pivrow (i))) ;
X}
X
Xdumpreport_ ()
X{
X fdumpstage (stdout) ;
X fdumpperm (stdout, "Row", & Pivrow (0), & Ipivrow (0)) ;
X fdumpperm (stdout, "Col", & Pivcol (0), & Ipivcol (0)) ;
X fdumpstats (stdout) ;
X}
X
Xfdumpreport_ ()
X{
X FILE *f ;
X char s [80], *p ;
X
X p = index (Name, ' ') ;
X if (p) *p = '\0' ;
X sprintf (s, "%s", Name) ;
X if (!Open_file (f, s, "a")) { printf ("Couldn't open %s\n", s) ; return ; }
X fdumpstage (f) ;
X fdumpperm (f, "Row", & Pivrow (0), & Ipivrow (0)) ;
X fdumpperm (f, "Col", & Pivcol (0), & Ipivcol (0)) ;
X fdumpstats (f) ;
X Close_file (f) ;
X}
X
Xfdumpstats (f)
XFILE *f ;
X{
X fprintf (f,
X "Initnz: %d L: %d U: %d L+U: %d Fillin: %d Error: %g\n",
X Nonzeros, Lnz, Unz, Lnz+Unz, Lnz+Unz-Nonzeros, Error) ;
X}
X
Xfdumpperm (f, what, piv, ipiv)
XFILE *f ;
Xchar *what ;
Xint *piv, *ipiv ;
X{
X int err = FALSE, i ;
X
X fprintf (f, "\n%s:\n", what) ;
X for (i = 1 ; i <= Size ; i++)
X {
X fprintf (f, " %4d", piv [i]) ;
X if (i % 15 == 0) fprintf (f, "\n") ;
X if (ipiv [piv [i]] != i) err = TRUE ;
X }
X fprintf (f, "\n") ;
X if (err) fprintf (f, "ERROR in I%s\n", what) ;
X}
X
Xfdumpstage (f)
XFILE *f ;
X{
X int i, count, prcol = 0 ;
X
X fprintf (f, "\n%s %5d stages:\n", Name, NStages) ;
X
X for (i = 1 ; i <= NStages ; i++)
X {
X if (Stagesize (i) == Stagesize (i+1))
X {
X count = 0 ;
X while (Stagesize (i) == Stagesize (i+1))
X {
X count++ ;
X i++ ;
X }
X if (prcol >= 70)
X {
X prcol = 0 ;
X fprintf (f, "\n") ;
X }
X fprintf (f, " %4d:%-4d", Stagesize (i), count) ;
X prcol += 10 ;
X }
X else
X {
X count = 0 ;
X if (prcol >= 75)
X {
X prcol = 0 ;
X printf ("\n") ;
X }
X fprintf (" %4d", Stagesize (i)) ;
X prcol += 5 ;
X }
X }
X}
SHAR_EOF
true || echo 'restore of dump.c failed'
fi
# ============= mark.c ==============
if test -f 'mark.c' -a X"$1" != X"-c"; then
echo 'x - skipping mark.c (File already exists)'
else
echo 'x - extracting mark.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'mark.c' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include "mattool.h"
X
Xmulti_proc ()
X{
X if (status == None || status == Done) load_proc () ;
X if (status == None || status == Done) return ;
X status = Running ;
X MFactor = dget_value (factor_item, 2.0) ;
X Sdense = dget_value (dense_item, -1.0) ;
X Idense = (int) Sdense ;
X printf ("MFactor: %f\n", MFactor) ;
X
X if (mstep == 2) m2 () ;
X if (mstep == 3) m3 () ;
X if (mstep == 4) m4 () ;
X while (Stage < Size) { m1 () ; m2 () ; m3 () ; m4 () ; }
X wrap_up () ;
X}
X
Xmstep_proc ()
X{
X if (status == None || status == Done) load_proc () ;
X if (status == None || status == Done) return ;
X status = Running ;
X MFactor = dget_value (factor_item, 2.0) ;
X Sdense = dget_value (dense_item, -1.0) ;
X Idense = (int) Sdense ;
X
X if (Stage >= Size) wrap_up () ;
X else if (mstep == 1) m1 () ;
X else if (mstep == 2) m2 () ;
X else if (mstep == 3) m3 () ;
X else if (mstep == 4) m4 () ;
X}
X
Xwrap_up ()
X{
X int row, col ;
X
X display_ALU () ;
X if (Printflag == 1) dumpmatrix_ () ;
X Lnz = Densesize * (Densesize-1) / 2 ;
X for (col = 1 ; col <= Sparsesize ; col++) Lnz += Llen (col) ;
X Unz = Densesize * (Densesize+1) / 2 ;
X for (row = 1 ; row <= Sparsesize ; row++) Unz += Ulen (row) ;
X mdone_ () ;
X/* fdumpreport_ () ; */
X fdumpstats (stdout) ;
X printf ("Sparsesize: %d Densesize: %d Idense: %d\n",
X Sparsesize, Densesize, Idense) ;
X mstep = 5 ;
X status = Done ;
X}
X
X
X/*----------------------------------------------------------------------------*/
X
Xm1 () /* sort the matrix according to the Markowitz data structures */
X{
X int i ;
X
X if (Densep != 0) return ;
X m1_ () ;
X if (Printflag == 1)
X {
X printf ("permutations:\n") ;
X for (i = 1 ; i <= Size ; i++)
X printf ("%3d %3d %3d %3d %3d %3d %3d\n", i,
X Pivrow (i), Ipivrow (Pivrow (i)), RowLen (Pivrow (i)),
X Pivcol (i), Ipivcol (Pivcol (i)), ColLen (Pivcol (i))) ;
X }
X if (display) display_ALU () ;
X mstep = 2 ;
X}
X
Xm2 () /* select a set of compatible pivots and highlight */
X{
X if (Densep != 0) return ;
X m2a_ () ;
X if (Densep != 0) return ;
X m2_ () ;
X if (display) display_pivots () ;
X mstep = 3 ;
X}
X
Xm3 () /* permute the pivots to the diagonal */
X{
X if (Densep != 0) return ;
X m3_ () ;
X if (display) display_ALU () ;
X printf ("Stage: %3d NPivots: %3d\n ", Stage, NPivots) ;
X for (i = 1 ; i <= NPivots ; i++) printf (" %3d", Cost (i)) ;
X printf ("\n") ;
X mstep = 4 ;
X}
X
Xm4 () /* perform the reductions and add fillin, update Markowitz data */
X{
X if (Densep != 0) return ;
X m4_ () ;
X if (Debug)
X {
X mywait ("A by cols") ;
X pw_batch_on (pw) ;
X display_clear () ;
X display_A () ;
X pw_batch_off (pw) ;
X mywait ("A by rows") ;
X pw_batch_on (pw) ;
X display_clear () ;
X display_A_rows () ;
X pw_batch_off (pw) ;
X mywait ("U") ;
X pw_batch_on (pw) ;
X display_clear () ;
X display_U () ;
X pw_batch_off (pw) ;
X mywait ("L") ;
X pw_batch_on (pw) ;
X display_clear () ;
X display_L () ;
X pw_batch_off (pw) ;
X mywait ("box") ;
X display_box () ;
X mywait ("update blocks") ;
X show_updates () ;
X mywait ("all") ;
X display_ALU () ;
X }
X else if (display) display_ALU () ;
X mstep = 1 ;
X}
X
Xmywait (c)
Xchar *c ;
X{
X int i ;
X printf ("Hit CR to see %s\n", c) ;
X i = getchar () ;
X}
X
Xcrwait_ ()
X{
X int i ;
X printf ("Hit CR to continue\n") ;
X i = getchar () ;
X}
SHAR_EOF
true || echo 'restore of mark.c failed'
fi
# ============= mattool.c ==============
if test -f 'mattool.c' -a X"$1" != X"-c"; then
echo 'x - skipping mattool.c (File already exists)'
else
echo 'x - extracting mattool.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'mattool.c' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include "mattool.h"
X
XFrame frame ;
XPanel control_panel ;
XCanvas canvas ;
XPanel_item dir_item, file_item, factor_item, dense_item ;
XMenu file_menu = NULL ;
Xint wr, wc, wsize, n, ratio, elsize ;
XPixwin *pw ;
X
Xstatic short icon_image[] = {
X#include "icon.h"
X} ;
X
XDEFINE_ICON_FROM_IMAGE (matrix_icon, icon_image) ;
X
XFILE *matfile ;
Xchar mat_title [80], key [10], type [10], path [256] ;
Xint i, reducs, step, ops, drops, init_nz ;
XDIR *cwd ;
Xchar pathname [MAXPATHLEN], filename [32] ;
Xchar files [Maxfiles][32], smsize [Maxfiles][8] ;
Xint msize [Maxfiles], entries = 0, status = Loaded, display = TRUE ;
Xint mstep = 1 ;
Xdouble *xm ;
X
Xglobal_data global_ ;
X
X/******************************************************************************/
X
Xcmain_ ()
X{
X static char *args [] = { "sparse_dynamics" } ;
X sparse_dynamics (1, args) ;
X}
X
Xsparse_dynamics (argc, argv)
Xint argc ;
Xchar **argv ;
X{
X int row, c, k ;
X Rect *fr, nfr ;
X
X Debug = FALSE ;
X Printflag = Debug ;
X strcpy (pathname, StartDir) ;
X
X /* start with an example 6-by-6 matrix */
X Size = 6 ;
X Nonzeros = 21 ;
X for (k = 1 ; k <= Nonzeros ; k++) A_value (k) = (double) k ;
X MemTail = Mmax+1 ;
X MemHead = Nonzeros+1 ;
X MemCTail = Cmax+1 ;
X MemCHead = Nonzeros+1 ;
X A_row ( 1) = 1 ; A_col ( 1) = 1 ;
X A_row ( 2) = 1 ; A_col ( 2) = 2 ;
X A_row ( 3) = 1 ; A_col ( 3) = 6 ;
X A_row ( 4) = 2 ; A_col ( 4) = 2 ;
X A_row ( 5) = 2 ; A_col ( 5) = 3 ;
X A_row ( 6) = 2 ; A_col ( 6) = 5 ;
X A_row ( 7) = 3 ; A_col ( 7) = 1 ;
X A_row ( 8) = 3 ; A_col ( 8) = 3 ;
X A_row ( 9) = 3 ; A_col ( 9) = 4 ;
X A_row (10) = 3 ; A_col (10) = 5 ;
X A_row (11) = 4 ; A_col (11) = 1 ;
X A_row (12) = 4 ; A_col (12) = 2 ;
X A_row (13) = 4 ; A_col (13) = 4 ;
X A_row (14) = 4 ; A_col (14) = 6 ;
X A_row (15) = 5 ; A_col (15) = 1 ;
X A_row (16) = 5 ; A_col (16) = 4 ;
X A_row (17) = 5 ; A_col (17) = 5 ;
X A_row (18) = 6 ; A_col (18) = 2 ;
X A_row (19) = 6 ; A_col (19) = 3 ;
X A_row (20) = 6 ; A_col (20) = 5 ;
X A_row (21) = 6 ; A_col (21) = 6 ;
X
X /* initialize canvas and control panel */
X
X frame = window_create (NULL, FRAME, FRAME_ARGS, argc, argv,
X FRAME_LABEL, "SparseDynamics Tool (with the D2 algorithm) Tim Davis \
X na.tdavis@na-net.ornl.gov", FRAME_ICON, &matrix_icon, 0) ;
X canvas = window_create (frame, CANVAS, WIN_COLUMNS, 120, WIN_ROWS, 57,
X CANVAS_MARGIN, 0, WIN_X, 0, 0) ;
X control_panel = window_create (frame, PANEL,
X WIN_RIGHT_OF, canvas, WIN_COLUMNS, 20, 0) ;
X
X row = -1 ;
X Button ("Step", mstep_proc, 0, ++row) ;
X Button ("Run", multi_proc, 0, ++row) ;
X Button ("Quit", quit_proc, 0, ++row) ;
X Toggle ("Display:", fast_proc, display, 0, ++row) ;
X Toggle ("Print:", print_proc, Printflag, 0, ++row) ;
X factor_item = Text ("Factor:", null_proc, 30, "2", 0, ++row) ;
X dense_item = Text ("Switch:", null_proc, 30, "-1", 0, ++row) ;
X dir_item = Text ("Dir: ", dir_proc, 30, pathname, 0, ++row) ;
X Button ("File:", file_proc, 0, ++row) ;
X
X file_item = panel_create_item (control_panel, PANEL_MESSAGE,
X PANEL_ITEM_X, ATTR_COL(9), PANEL_ITEM_Y, ATTR_ROW(row), 0) ;
X window_set (control_panel, PANEL_CARET_ITEM, dir_item, 0) ;
X dir_setup () ;
X
X window_fit (frame) ;
X fr = (Rect *) window_get (frame, FRAME_OPEN_RECT) ;
X nfr.r_left = 0 ;
X nfr.r_top = 0 ;
X nfr.r_width = fr -> r_width ;
X nfr.r_height = fr -> r_height ;
X window_set (frame, FRAME_OPEN_RECT, &nfr, 0) ;
X
X /* store and display original matrix */
X
X scale_display () ;
X display_original () ;
X storematrix_ () ;
X
X /* relinquish control to SunTools */
X
X window_main_loop (frame) ;
X exit (0) ;
X}
X
Xnull_proc () { }
SHAR_EOF
true || echo 'restore of mattool.c failed'
fi
# ============= fdefine.h ==============
if test -f 'fdefine.h' -a X"$1" != X"-c"; then
echo 'x - skipping fdefine.h (File already exists)'
else
echo 'x - extracting fdefine.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fdefine.h' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X/*----------------------------------------------------------------------------*/
X/* A' matrix during factorization: */
X
X#define Row_next(p) Im(p)
X#define Row_last(p) Xm(p) /* requires (int) or idint ( ) */
X#define Row_index(p,j) Im(p+j)
X#define Row_value(p,j) Xm(p+j)
X
X#define Address_of_Row_index(p,j) (p+j)
X
X#define Col_next(p) Cm(p)
X#define Col_index(p,j) Cm(p+j)
X#define Col_last(p) Cm(p+Bl1)
X
X#define Address_of_Col_index(p,j) (p+j)
X
X/*----------------------------------------------------------------------------*/
X/* L and U factors during and after factorization: */
X
X#define L_index(p,j) Im(p+j-1)
X#define L_value(p,j) Xm(p+j-1)
X
X#define Address_of_L_index(p,j) (p+j-1)
X
X#define L_next(p) Im(p)
X#define L_col(p) Xm(p) /* requires (int) or idint ( ) */
X
X#define U_index(p,j) Im(p+j-1)
X#define U_value(p,j) Xm(p+j-1)
X
X/*----------------------------------------------------------------------------*/
X/* original A matrix before factorization: */
X
X#define A_row(k) Cm(k)
X#define A_col(k) Im(k)
X#define A_value(k) Xm(k)
X
X/*----------------------------------------------------------------------------*/
X/* Rowia, Colia structure: */
X
X#define Row_ia(p,j) Im(p+j-1)
X#define Col_ia(p,j) Im(p+j-1)
SHAR_EOF
true || echo 'restore of fdefine.h failed'
fi
# ============= global.h ==============
if test -f 'global.h' -a X"$1" != X"-c"; then
echo 'x - skipping global.h (File already exists)'
else
echo 'x - extracting global.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'global.h' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
Xc constants:
Xc
Xc Nil the Nil pointer value (0)
Xc
Xc Mmax size of integer free space (im), size of real*8 free space (xm)
Xc Cmax size of integer free space for col. storage during factorization
Xc Nmax maximum matrix size
Xc
Xc Bl
Xc Rowbl
Xc Colbl
X
Xc IF YOU CHANGE THESE PARAMETETERS, BE SURE TO CHANGE CDEFINE.H !!!!!
X integer Mmax, Cmax, Nmax, Nil, Bl, Bl1, Rowbl, Colbl
X parameter (Mmax = 200000, Cmax=180000, Nmax = 823, Nil=0)
X parameter (Bl=32, Bl1=Bl+1, Rowbl=Bl+1, Colbl=Bl+2)
X
Xc-------------------------------------------------------------------------------
Xc Allocatable memory:
Xc
Xc Im integer free space
Xc Xm real free space
Xc
Xc MemHead index of first free space in im/xm
Xc MemTail index of highest used space in im/xm
Xc MemCHead index of first free space in cm
Xc MemCTail index of highest used space in cm
X
X integer Im (Mmax), MemHead, MemTail, MemLock
X integer Cm (Cmax), MemCHead, MemCTail, MemCLock
X real*8 Xm (Mmax)
X
Xc-------------------------------------------------------------------------------
Xc Data structure for input matrix A:
Xc
Xc Size size of A
Xc Nonzeros number of nonzeros in A
X
X integer Size, Nonzeros
X character Name*80
X
Xc-------------------------------------------------------------------------------
Xc Data structure for A':
Xc
Xc RowLen number of nonzero elements in rows of A'
Xc RowHead pointers (in Im) to rows of A'
Xc forms link list of Bl-element blocks with indices and values
Xc
Xc RowHead record for row i (where first block p = RowHead (i)):
Xc integer next | Im (p), pointer to next block
Xc | last block holds Nil
Xc integer last | Xm (p), pointer to prev. block
Xc | first block holds -i
Xc integer index (1..Bl) | Im (p+1 .. p+Bl)
Xc real*8 value (1..Bl) | Xm (p+1 .. p+Bl)
Xc | next block starts at Im (p + Rowbl)
Xc
Xc ColLen number of nonzero elements in columns of A'
Xc ColHead pointers (in Cm) to columns of A'
Xc forms link list of Bl-element blocks with indices only
Xc
Xc ColHead record for col i (where first block p = ColHead (i)):
Xc integer next | Cm (p), pointer to next block
Xc | last block holds Nil
Xc integer index (1..Bl) | Cm (p+1 .. p+Bl)
Xc integer last | Cm (p+Bl+1), pointer to prev. block
Xc | first block holds -i
Xc | next block starts at Cm (p + Colbl)
X
X integer RowLen (Nmax), RowHead (Nmax)
X integer ColLen (Nmax), ColHead (Nmax)
X
Xc-------------------------------------------------------------------------------
Xc Data structure for L and U:
Xc
Xc Lp pointers (in Im) to columns of L
Xc Llen length of each column of L
Xc
Xc Up pointer (in Im) to rows of U
Xc Ulen length of each row of U
Xc
Xc Structure for a single column/row of L/U (where p = Lp (i) or Up (i)
Xc and len = Llen (i) or Ulen (i)):
Xc
Xc integer index (1..len) | Im (p .. p+len-1)
Xc real*8 value (1..len) | Xm (p .. p+len-1)
Xc | next block starts at Im (p + len)
X
X integer Lp (Nmax), Up (Nmax), Llen (Nmax), Ulen (Nmax)
X
Xc-------------------------------------------------------------------------------
Xc Permutation matrices:
Xc
Xc Pivrow Pivrow (new) gives old row index
Xc Pivcol Pivcol (new) gives old column index
Xc Ipivrow Ipivrow (old) gives new row index
Xc Ipivrow Ipivcol (old) gives new column index
X
X integer Pivrow (Nmax), Ipivrow (Nmax)
X integer Pivcol (Nmax), Ipivcol (Nmax)
X
Xc-------------------------------------------------------------------------------
Xc Other:
Xc
Xc Iw1 - Iw4 integer work spaces
Xc W real work space
Xc Wi integer work space
Xc Rhs right-hand side of Ax=b
Xc Xt true solution
Xc Error relative infinite error: max (x - x~) / max (x)
Xc Prows Prows (1 .. NPivots) gives pivot rows at each stage
Xc Pcols Pcols (1 .. NPivots) gives pivot columns at each stage
Xc NStages number of parallel stages so far
Xc Stage number of pivots processed so far
Xc Printflag 1 if printing, 0 if not
Xc Stagesize number of pivots in each parallel step
Xc Pfirst Stage + 1
Xc Plast Stage + NPivots
Xc Redrow reduce rows Redrow (1..NRedrow) in A'
Xc Redcol reduce cols Redcol (1..NRedcol) in A'
Xc Lpointer head of link list for each row to reduce (equiv. to Pcols)
Xc Debug 1 if debuging, 0 if not
X
X integer Prows (Nmax), Pcols (Nmax), Stagesize (Nmax)
X integer Iw1 (Nmax+2), Iw2 (Nmax+2), Iw3 (Nmax+2), Iw4 (Nmax+2)
X integer Stage, NStages, Printflag, NPivots, Pfirst, Plast
X integer Redrow (Nmax), Redcol (Nmax), NRedrow, NRedcol, Debug
X integer Lpointer (Nmax), Wi (Nmax+2)
X equivalence (Lpointer, Pcols)
X real*8 Rhs (Nmax), Xt (Nmax), W (Nmax+2), Error
X integer Lnz, Unz
X real*8 Sdense
X integer Idense, Densesize, Densep, Sparsesize
X
Xc-------------------------------------------------------------------------------
Xc Markowitz selection data structure:
Xc
Xc MRowhead (nz) pointer to list of rows with nz nonzeros
Xc MRownext pointer to next row
Xc MRowlast pointer to last row
Xc
Xc MColhead (nz) pointer to list of columns with nz nonzeros
Xc MColnext pointer to next columns
Xc MCollast pointer to last columns
Xc
Xc MRowmin minimum number of nonzeros in any row
Xc MColmin minimum number of nonzeros in any column
Xc
Xc MRowflag 1 if row is selected, 0 if not
Xc MColflag 1 if col is selected, 0 if not
Xc
Xc MRowmax abs. maximum value in row (-1 if not computed)
Xc MFactor parallelism factor for Markowitz selection
X
X integer MRowhead (Nmax), MRownext (Nmax), MRowlast (Nmax)
X integer MColhead (Nmax), MColnext (Nmax), MCollast (Nmax)
X integer MRowmin, MColmin, MRowflag (Nmax), MColflag (Nmax)
X real*8 MRowmax (Nmax), MFactor
X
Xc-----------------------------------------------------------------------
X
X common /global/
Xc real*8
X 1 Xm, Rhs, Xt, W, MRowmax, MFactor, Error, Sdense,
Xc integer
X 1 Im, MemHead, MemTail, MemLock,
X 1 Cm, MemCHead, MemCTail, MemCLock,
X 1 Size, Nonzeros,
X 1 RowLen, RowHead,
X 1 ColLen, ColHead,
X 1 Lp, Up, Llen, Ulen,
X 1 Pivrow, Ipivrow, Pivcol, Ipivcol,
X 1 MRowhead, MRownext, MRowlast,
X 1 MColhead, MColnext, MCollast,
X 1 MRowmin, MColmin, MRowflag, MColflag,
X 1 Prows, Pcols, Stagesize,
X 1 Iw1, Iw2, Iw3, Iw4, Wi,
X 1 Stage, NStages, Printflag, NPivots,
X 1 Pfirst, Plast, Debug, Idense, Densesize, Densep, Sparsesize,
X 1 Redrow, Redcol, NRedrow, NRedcol, Lnz, Unz,
Xc character
X 1 Name
SHAR_EOF
true || echo 'restore of global.h failed'
fi
# ============= m1.F ==============
if test -f 'm1.F' -a X"$1" != X"-c"; then
echo 'x - skipping m1.F (File already exists)'
else
echo 'x - extracting m1.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'm1.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
Xc sort the matrix according to the Markowitz data structures
Xc (this step is for viewing purposes only)
X
X subroutine m1
X
X include 'global.h'
X integer nz, row, col, i
X
X if (Printflag .eq. 1) print *, 'm1', Stage, Size
X
Xc sort matrix according to markowitz selection data structure
Xc (rows Stage+1 to Size and cols Stage+1 to Size)
X
X if (MRowmin .eq. 0 .or. MColmin .eq. 0) return
X
X nz = MRowmin
X row = MRowhead (nz)
X do 20 i = Stage+1, Size-1
X Ipivrow (row) = i
X Pivrow (i) = row
X row = MRownext (row)
X if (row .eq. Nil) then
X10 nz = nz + 1
X if (nz .gt. Size) then
X print *, 'ROW NZ > SIZE!'
X call crwait
X return
X endif
X row = MRowhead (nz)
X if (row .eq. Nil) goto 10
X endif
X20 continue
X Ipivrow (row) = Size
X Pivrow (Size) = row
X
X nz = MColmin
X col = MColhead (nz)
X do 40 i = Stage+1, Size-1
X Ipivcol (col) = i
X Pivcol (i) = col
X col = MColnext (col)
X if (col .eq. Nil) then
X30 nz = nz + 1
X if (nz .gt. Size) then
X print *, 'COL NZ > SIZE!'
X call crwait
X return
X endif
X col = MColhead (nz)
X if (col .eq. Nil) goto 30
X endif
X40 continue
X Ipivcol (col) = Size
X Pivcol (Size) = col
X
X return
X end
SHAR_EOF
true || echo 'restore of m1.F failed'
fi
# ============= m2.F ==============
if test -f 'm2.F' -a X"$1" != X"-c"; then
echo 'x - skipping m2.F (File already exists)'
else
echo 'x - extracting m2.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'm2.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
Xc select a set of compatible pivots
X
X subroutine m2
X
X include 'global.h'
X integer prow, pcol, cost, mincost, maxcost, minnz, pp, kk, nz
X integer row, col, prownz, pcolnz, large, ik, ip, p, k, maxnz
X integer i, j, oldrow, oldcol
X logical pr, db
X real*8 findrowmax, amax, au
X
X pr = Printflag .eq. 1
X db = Debug .eq. 1
X NPivots = 0
X NRedrow = 0
X NRedcol = 0
X large = size*size + 1
X
X if (pr) then
X print *, 'm2'
X print *,'------------------------------------------------------'
X endif
X
X if (MFactor .lt. 0) then
X call markowitz (mincost, prow, pcol)
X if (mincost .lt. 0) return
X call pivotfound (prow, pcol, row, col, 3)
X if (pr) print *, 'SINGLE pivot:', mincost, prow, pcol
X return
X endif
X
Xc find all col pivots with cost = 0
X
X if (MColmin .eq. 1) then
X col = MColhead (1)
X30 if (col .eq. Nil) goto 40
X row = Col_index (ColHead (col), 1)
X if (MRowflag (row) .eq. 0) then
X prow = row
X pcol = col
X call pivotfound (prow, pcol, row, col, -2)
X if (pr) print *, 'col zero pivot:', prow, pcol
X else
X if (pr) print *, 'col zero pivot:', row, col,
X 1 ' incompatible'
X col = MColnext (col)
X endif
X goto 30
X endif
X40 continue
X
Xc find all row pivots with cost = 0
X
X if (MRowmin .eq. 1) then
X row = MRowhead (1)
X10 if (row .eq. Nil) goto 50
X col = Row_index (RowHead (row), 1)
X if (MColflag (col) .eq. 0) then
X prow = row
X pcol = col
X call pivotfound (prow, pcol, row, col, -1)
X if (pr) print *, 'row zero pivot:', prow, pcol
X else
X if (pr) print *, 'row zero pivot:', row, col,
X 1 ' incompatible'
X row = MRownext (row)
X endif
X goto 10
X endif
X50 if ((pr).and.NPivots.gt.0) print *,'zero cost pivots: ',NPivots
X
Xc find a single pivot with lowest nonzero Markowitz cost
X call markowitz (mincost, prow, pcol)
X if (mincost .lt. 0) return
X call pivotfound (prow, pcol, row, col, 3)
X if (pr) print *, 'min cost pivot:', mincost, prow, pcol
X
Xc find a set of compatible pivots with nonzero cost
X
X maxcost = MFactor * mincost
X minnz = max (2, min (MRowmin, MColmin))
X maxnz = int (sqrt (float (maxcost))) + 1
X if (pr) then
X print *, 'mincost: ', mincost, ' maxcost:', maxcost
X print *, 'minnz:',minnz,' maxnz:',maxnz
X endif
X
X do 400 nz = minnz, maxnz
X
Xc scan rows for the lowest cost pivot in the rows
X
X row = MRowhead (nz)
X200 if (row .eq. Nil) goto 230
X amax = MRowmax (row)
X au = amax * 0.01d0
X pcolnz = large
X pcol = 0
X k = 1
X p = RowHead (row)
X do 220 i = 1, nz
X col = Row_index (p, k)
X if (MColflag (col) .eq. 1) then
X if (db) then
X print *, 'incompatible row pivot:', row, col
X endif
X elseif (ColLen (col) .lt. pcolnz) then
X if (amax .lt. 0.0d0) then
X amax = findrowmax (row, col, ik, ip)
X MRowmax (row) = amax
X au = amax * 0.01d0
X endif
X if (abs (Row_value (p, k)) .gt. au) then
X pcol = col
X pcolnz = ColLen (col)
X endif
X endif
X call nextrow (p, k)
X220 continue
X prow = row
X cost = (nz - 1) * (pcolnz - 1)
X if (cost .le. maxcost .and. pcol .ne. 0) then
X if (pr) print *, 'row pivot', cost, prow, pcol
X oldrow = row
X call pivotfound (prow, pcol, row, col, 4)
X if (oldrow .eq. row) then
X print *, 'WARNING: oldrow .eq. row', row
X call crwait
X endif
X else
X row = MRownext (row)
X endif
X goto 200
X230 continue
X
Xc scan cols for the lowest cost pivot in the col
X
X col = MColhead (nz)
X if ((db) .and. col .ne. Nil) then
X print *,'::::::::::::::::::::::::: start testing colums'
X endif
X
X300 if (col .eq. Nil) goto 350
X if (db) then
X print *, '------------------------------------------'
X print *, 'testing column ', col, ' nz=', nz
X call dumpcol (col)
X print *, 'checking each element in column:'
X endif
X prownz = large
X prow = 0
X k = 1
X p = ColHead (col)
X do 340 i = 1, nz
X row = Col_index (p, k)
X if (MRowflag (row) .eq. 1) then
X if (db) then
X print *, 'incompatible col pivot:', row, col, p, k
X endif
X elseif (RowLen (row) .lt. prownz) then
X if (db) print *, 'testing col piv: ', row, col, p, k
X amax = MRowmax (row)
X if (amax .lt. 0.0d0) then
X amax = findrowmax (row, col, ik, ip)
X MRowmax (row) = amax
X else
X kk = 1
X pp = RowHead (row)
X do 320 j = 1, RowLen (row)
X if (col .eq. Row_index (pp, kk)) then
X ik = kk
X ip = pp
X goto 330
X endif
X call nextrow (pp, kk)
X320 continue
X endif
X330 au = amax * 0.01d0
X
X if (db) then
X print *, 'candidate: ', row
X print *, 'index: ', Row_index (ip,ik)
X print *, 'value: ', Row_value (ip,ik)
X print *, 'ip, ik:', ip,ik, ' row nz:',RowLen(row)
X print *, 'amax:', amax, ' au:', au
X endif
X
X if (abs (Row_value (ip, ik)) .gt. au) then
X prow = row
X prownz = RowLen (row)
X if (db) then
X print *, 'THIS ONE IS OK!'
X print *, 'maxcost:', maxcost
X print *, 'cost of this one:',(nz-1)*(prownz-1)
X endif
X endif
X endif
X call nextcol (p, k)
X340 continue
X pcol = col
X cost = (nz - 1) * (prownz - 1)
X if (cost .le. maxcost .and. prow .ne. 0) then
X if (pr) print *, 'col pivot cost:', cost, prow, pcol
X oldcol = col
X call pivotfound (prow, pcol, row, col, 5)
X if (oldcol .eq. col) then
X print *, 'WARNING: oldcol .eq. col', col
X call crwait
X endif
X else
X col = MColnext (col)
X endif
X goto 300
X350 continue
X
X if (pr) then
X print *,':::::::::::::::::::::::done testing colums'
X endif
X
X400 continue
X
X return
X end
X
X real*8 function findrowmax (row, col, ik, ip)
X include 'global.h'
X integer row, col, p, k, i, len, ik, ip
X real*8 amax
X
X amax = 0.0d0
X if (Debug .eq. 1) print *, 'findrowmax:', row, col
X len = RowLen (row)
X p = RowHead (row)
X k = 1
X do 10 i = 1, len
X amax = max (amax, abs (Row_value (p, k)))
X if (Debug.eq.1)print *,p,k,Row_index(p,k),Row_value(p,k),amax
X if (Row_index (p, k) .eq. col) then
X ik = k
X ip = p
X endif
X call nextrow (p, k)
X 10 continue
X findrowmax = amax
X return
X end
X
X subroutine nextrow (p, k)
X include 'global.h'
X integer p, k
X k = k + 1
X if (k .gt. Bl) then
X k = 1
X p = Row_next (p)
X endif
X return
X end
X
X subroutine nextcol (p, k)
X include 'global.h'
X integer p, k
X k = k + 1
X if (k .gt. Bl) then
X k = 1
X p = Col_next (p)
X endif
X return
X end
SHAR_EOF
true || echo 'restore of m2.F failed'
fi
# ============= m2a.F ==============
if test -f 'm2a.F' -a X"$1" != X"-c"; then
echo 'x - skipping m2a.F (File already exists)'
else
echo 'x - extracting m2a.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'm2a.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
Xc check for swith to dense
X
X subroutine m2a
X include 'global.h'
X integer need
X
X if (Idense .lt. 0) return
X if (NStages + 1 .lt. Idense) return
X
X print *, 'switching to dense'
X Sparsesize = Stage
X Densesize = Size - Stage
X print *, 'sparse size:', Sparsesize, ' dense size:', Densesize
X
Xc move sparse A' into dense storage
X
X need = Densesize**2
X Densep = MemHead
X MemHead = MemHead + need
X if (MemTail .lt. MemHead) then
X print *, 'Out of Memory for dense matrix!', MemHead, MemTail
X return
X endif
X
Xc factor dense A' into L and U
X
X call denselu (Xm (Densep), Densesize, Sparsesize)
X Stage = Size
X return
X end
X
X subroutine denselu (a, n, s)
X include 'global.h'
X real*8 a (n, n), swap, alpha, amax
X integer row, col, oldrow, p, len, j, k, n, kl, s
X integer pp, kk
X
Xc call dumpmatrix
X
X do 20 row = 1, n
X do 10 col = 1, n
X10 a (row, col) = 0.0d0
X oldrow = Pivrow (row + s)
X p = RowHead (oldrow)
X len = RowLen (oldrow)
X do 20 j = 1, len, Bl
X kl = min (Bl, len - j + 1)
X do 30 k = 1, kl
X col = Ipivcol (Row_index (p, k)) - s
X30 a (row, col) = Row_value (p, k)
X20 p = Row_next (p)
X
Xc print *, 'initial matrix:'
Xc call dumpdense (a, n, s, Pivrow, Ipivrow)
X
X do 90 k = 1, n-1
X
Xc find the pivot row
X
X amax = abs (a (k, k))
X p = k
X do 60 row = k+1, n
X if (amax .lt. abs (a (row, k))) then
X amax = abs (a (row, k))
X p = row
X endif
X60 continue
X
Xc swap the pivot row with row k
X
X if (k .ne. p) then
X do 70 col = 1, n
X swap = a (p, col)
X a (p, col) = a (k, col)
X70 a (k, col) = swap
X
X pp = Pivrow (p + s)
X kk = Pivrow (k + s)
X Pivrow (p + s) = kk
X Pivrow (k + s) = pp
X Ipivrow (pp) = k + s
X Ipivrow (kk) = p + s
X endif
X
Xc print *, 'pivot', k
Xc call dumpdense (a, n, s, Pivrow, Ipivrow)
X
X if (a (k, k) .eq. 0.0d0) then
X print *, 'ERROR: alpha is undefined!', k
X return
X endif
X
Xc reduce the matrix using the pivot row
X
X do 80 row = k+1, n
X alpha = a (row, k) / a (k, k)
X a (row, k) = alpha
X do 80 col = k+1, n
X80 a (row, col) = a (row, col) - alpha * a (k, col)
X
Xc print *, 'reduced', k
Xc call dumpdense (a, n, s, Pivrow, Ipivrow)
X90 continue
X return
X end
X
X subroutine dumpdense (a, n, s, p, ip)
X include 'global.h'
X real*8 a (n,n)
X integer n, row, col, p (Size), ip (Size), s
X do 22 row = 1, n
X print *, ' '
X22 print 11, row, p (row+s), ip (p (row+s)), (a (row, col), col = 1, n)
X11 format (' ', i3, i3, i3,': ', /, (9g12.3))
X print *, ' '
X return
X end
SHAR_EOF
true || echo 'restore of m2a.F failed'
fi
# ============= m3.F ==============
if test -f 'm3.F' -a X"$1" != X"-c"; then
echo 'x - skipping m3.F (File already exists)'
else
echo 'x - extracting m3.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'm3.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
Xc permute the pivots to the diagonal
X
X subroutine m3
X
X include 'global.h'
X integer pcol, prow, p, pp, i, ip, m, col, row
X
X integer j, temp
X
X if (Printflag .eq. 1) print *, 'm3'
X
X Pfirst = Stage + 1
X Plast = Stage + NPivots
X
Xc------------------------------------------------------------------------------
Xc sort pivots according to original ordering (sort by Prows, then Pcols
Xc follows)
X
X do 5 i = 1, NPivots
X
X do 4 j = i, NPivots
X if (Prows (i) .gt. Prows (j)) then
X temp = Prows (i)
X Prows (i) = Prows (j)
X Prows (j) = temp
X
X temp = Pcols (i)
X Pcols (i) = Pcols (j)
X Pcols (j) = temp
X endif
X4 continue
X5 continue
X
Xc------------------------------------------------------------------------------
X
X
X do 10 i = 1, NPivots
X
X pcol = Pcols (i)
X prow = Prows (i)
X p = Stage + i
X
X pp = Pivcol (p)
X ip = Ipivcol (pcol)
X Pivcol (ip) = pp
X Pivcol (p) = pcol
X Ipivcol (pp) = ip
X Ipivcol (pcol) = p
X
X pp = Pivrow (p)
X ip = Ipivrow (prow)
X Pivrow (ip) = pp
X Pivrow (p) = prow
X Ipivrow (pp) = ip
X 10 Ipivrow (prow) = p
X
Xc at this point, Prows and Pcols are no longer needed
X
X NStages = NStages + 1
X Stagesize (NStages) = NPivots
X
Xc-------
Xc permute rows and cols of A' that will be modified to upper left part of A'
Xc (this is for viewing purposes only)
X
X do 120 i = Plast + 1, Size
X 120 Pcols (i) = Pivcol (i)
X
X do 130 i = 1, NRedcol
X m = Plast + i
X Pivcol (m) = Redcol (i)
X 130 Ipivcol (Pivcol (m)) = m
X
X m = Plast + NRedcol
X do 140 i = Plast + 1, Size
X col = Pcols (i)
X if (MColflag (col) .eq. 0) then
X m = m + 1
X Pivcol (m) = col
X Ipivcol (Pivcol (m)) = m
X endif
X 140 continue
X
X if (m .ne. Size) then
X print *, 'ERROR: m .ne. Size', m, Size, NRedcol, NRedrow
X print *, ' in col'
X call crwait
X return
X endif
X
X do 220 i = Plast + 1, Size
X 220 Prows (i) = Pivrow (i)
X
X do 30 i = 1, NRedrow
X m = Plast + i
X Pivrow (m) = Redrow (i)
X 30 Ipivrow (Pivrow (m)) = m
X
X m = Plast + NRedrow
X do 40 i = Plast + 1, Size
X row = Prows (i)
X if (MRowflag (row) .eq. 0) then
X m = m + 1
X Pivrow (m) = row
X Ipivrow (Pivrow (m)) = m
X endif
X 40 continue
X
X if (m .ne. Size) then
X print *, 'ERROR: m .ne. Size', m, Size, NRedcol, NRedrow
X print *, ' in row'
X call crwait
X return
X endif
X
X if (Printflag .eq. 1) then
X print *, 'NRedcol: ', NRedcol
X print *, 'NRedrow: ', NRedrow
X endif
Xc-------
X return
X end
SHAR_EOF
true || echo 'restore of m3.F failed'
fi
# ============= m4.F ==============
if test -f 'm4.F' -a X"$1" != X"-c"; then
echo 'x - skipping m4.F (File already exists)'
else
echo 'x - extracting m4.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'm4.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
X subroutine m4
X
X include 'global.h'
X integer pcol, prow, p, i, len, k, j, ap, row, plen, pl
X integer pu, piv, kl, lastp, nextp, Fi (Nmax), nfill, nz
X integer need, needblocks, npiv, npivs, sp, sk, col, b, l
X integer pivlink, nextlink, memsize, lastblock, nextblock
X real*8 alpha, zz, Fx (Nmax)
X equivalence (Fi, Iw1)
X logical pr, db
X
X pr = Printflag .eq. 1
X db = Debug .eq. 1
X
X if (pr) print *, 'm4'
X if (db) call dumpmatrix
X Stage = Stage + NPivots
X
Xc pivots have been chosen (Pivcol (Pfirst..Plast) and Pivrow (Pfirst..Plast))
X
Xc allocate space for new columns of L
X
X do 10 i = Pfirst, Plast
X pcol = Pivcol (i)
X len = ColLen (pcol) - 1
X p = MemHead
X MemHead = MemHead + len
X Lp (i) = p
X 10 Llen (i) = len
X
Xc allocate space for new rows of U
X
X do 20 i = Pfirst, Plast
X prow = Pivrow (i)
X len = RowLen (prow)
X p = MemHead
X MemHead = MemHead + len
X Up (i) = p
X 20 Ulen (i) = len
X
X if (MemTail .lt. MemHead) then
X print *, 'Out of Im/Xm memory for LU!', MemHead, MemTail
X return
X endif
X
Xc store rows Pivrow (Pfirst .. Plast) of A' into new rows of U
Xc (putting pivot element first)
X
X do 40 i = Pfirst, Plast
X prow = Pivrow (i)
X pcol = Pivcol (i)
X ap = RowHead (prow)
X k = 1
X p = Up (i)
X len = Ulen (i)
X do 40 j = 1, len
X if (Row_index (ap, k) .eq. pcol) then
X U_index (p, j) = U_index (p, 1)
X U_value (p, j) = U_value (p, 1)
X U_index (p, 1) = Row_index (ap, k)
X U_value (p, 1) = Row_value (ap, k)
X else
X U_index (p, j) = Row_index (ap, k)
X U_value (p, j) = Row_value (ap, k)
X endif
X call nextrow (ap, k)
X 40 continue
X
Xc store columns Pivcol (Pfirst .. Plast) of A' into new columns of L
Xc (removing pivot element)
Xc also set up pointer chain through L for row to reduce
X
X do 50 i = 1, NRedrow
X 50 Lpointer (Redrow (i)) = Nil
X
X do 60 i = Pfirst, Plast
X prow = Pivrow (i)
X pcol = Pivcol (i)
X ap = ColHead (pcol)
X k = 1
X p = Lp (i)
X len = Llen (i)
X do 60 j = 1, len
X if (Col_index (ap, k) .eq. prow) then
X call nextcol (ap, k)
X endif
X row = Col_index (ap, k)
X L_index (p, j) = Lpointer (row)
X Lpointer (row) = Address_of_L_index (p, j)
X L_value (p, j) = i
X call nextcol (ap, k)
X 60 continue
X
X if (db) then
X print *, 'After LU factors have been formed:'
X call dumpu
X call dumpllist
X call dumpmatrix
X endif
X
Xc-------------------------------------------------------------------------------
Xc for each row to modify: (parallel loop)
X
X do 100 i = 1, NRedrow
X
X row = Redrow (i)
X p = RowHead (row)
X len = RowLen (row)
X MRowmax (row) = -1.0d0
X
Xc scatter row R to be reduced into W
X
X if (db) then
X print *,'-----------------------------------------------'
X print *,'Reducing row ',row,' at i=',i,' p=',p,' len=',len
X print *,'Ipivrow (row)=', Ipivrow (row)
X print *, 'scatter R into W'
X endif
X
X pivlink = Nil
X npivs = 0
X kl = 0
X do 110 j = 1, len, Bl
X kl = min (Bl, len - j + 1)
X do 120 k = 1, kl
X col = Row_index (p, k)
X W (col) = Row_value (p, k)
X Wi (col) = 1
X if (Ipivcol (col) .le. Plast) then
X Row_index (p, k) = pivlink
X pivlink = Address_of_Row_index (p, k)
X npivs = npivs + 1
X if (db) print 95, j, p, k, col, Ipivcol (col),
X 1 W (col), pivlink, npivs
X else
X if (db) print 97, j, p, k, col, Ipivcol (col),W(col)
X endif
X 120 continue
X lastp = p
X 110 p = Row_next (p)
X
Xc eliminate elements in pivot columns from R
X
X if (db) print *, 'eliminate pivots from R'
X p = lastp
X k = kl
X 115 nextlink = Im (pivlink)
X Im (pivlink) = Row_index (p, k)
X Xm (pivlink) = Row_value (p, k)
X if (db) print *, 'pivlink=', pivlink, ' p,k', p,k, ' i,x',
X 1 Im (pivlink), Xm (pivlink)
X pivlink = nextlink
X if (pivlink .eq. Nil) goto 125
X k = k - 1
X if (k .le. 0) then
X k = Bl
X p = idint (Row_last (p))
X endif
X goto 115
X 125 len = len - npivs
X RowLen (row) = len
X if (db) call dumprow (row)
X
Xc scatter/add each pivot row P1..Pk to W
Xc store row and alpha into L data structure
X
X p = Lpointer (row)
X npiv = 0
X
X 140 npiv = npiv + 1
X piv = idint (L_col (p))
X nextp = L_next (p)
X pu = Up (piv)
X plen = Ulen (piv)
X if (db) then
X print *, ' '
X print *, 'add piv=',piv,' pu=',pu,' plen=',plen,' p=',p
X endif
X if (Pivcol (piv) .ne. U_index (pu, 1)) then
X print *, 'ERROR Piv', Pivcol (piv), U_index (pu, 1)
X return
X endif
X alpha = W (Pivcol (piv)) / U_value (pu, 1)
X Fi (npiv) = piv
X Fx (npiv) = alpha
X if (db) then
X print *, 'alpha=',alpha, ' Pivcol(piv)=',Pivcol(piv)
X print *, 'W (Pivcol (piv)) =', W (Pivcol (piv))
X print *, 'U_value (pu,1) =',U_value(pu,1),' npiv=',npiv
X endif
X W (Pivcol (piv)) = 0.0d0
X Wi (Pivcol (piv)) = 0
X if (alpha .eq. 0.0d0) then
X if (db) then
X print *, 'WARNING: alpha is zero'
Xc call crwait
X endif
X endif
X if (U_value (pu, 1) .eq. 0.0d0) then
X print *, 'ERROR: alpha is undefined!'
X return
X endif
X do 130 k = 2, plen
X col = U_index (pu, k)
X zz = W (col) - alpha * U_value (pu, k)
X if (db) then
X print 99, col, Ipivcol (col), zz, W (col), Wi (col),
X 1 alpha, U_value(pu,k)
X endif
X Wi (col) = 1
X 130 W (col) = zz
X Im (p) = row
X Xm (p) = alpha
X p = nextp
X if (p .ne. Nil) goto 140
X
X if (db) print *, 'number of pivot rows added to R: ', npiv
X if (npiv .ne. npivs) then
X print *, 'ERROR npiv .ne. npivs', npiv, npivs
X return
X endif
X
Xc gather nonzeros in R back from W using Rowindex, zeroing W as we go
X
X if (db) print *, 'Gather nonzeros in R, zeroing W as we go'
X p = RowHead (row)
X kl = 0
X lastp = p
X
X do 150 j = 1, len, Bl
X kl = min (Bl, len - j + 1)
X do 160 k = 1, kl
X col = Row_index (p, k)
X if (db) print 97, j, p, k, col, Ipivcol (col), W (col)
X Row_value (p, k) = W (col)
X if (W (col) .eq. 0.0d0) then
X if (db) then
X print *, 'WARNING: zero as update'
Xc call crwait
X endif
X endif
X W (col) = 0.0d0
X Wi (col) = 0
X 160 continue
X lastp = p
X 150 p = Row_next (p)
X
X sp = lastp
X sk = kl
X if (db) print *, 'nonzeros so far ', len, ' sp, sk', sp, sk
X
Xc gather nonzeros still in W into Fi and Fx
X
X nfill = 0
X
X do 170 npiv = 1, npivs
X piv = Fi (npiv)
X alpha = Fx (npiv)
X pu = Up (piv)
X plen = Ulen (piv)
X if (db) then
X print *, ' '
X print *, 'fillin piv=',piv,' pu=',pu
X 1 , ' plen=',plen,' p=',p
X endif
X
X do 180 k = 2, plen
X col = U_index (pu, k)
X if (Wi (col) .ne. 0) then
X nfill = nfill + 1
X Fi (npivs + nfill) = col
X Fx (npivs + nfill) = W (col)
X if (W (col).eq.0.0d0) then
X if (db) then
X print *,'WARNING: zero as fillin'
Xc call crwait
X endif
X endif
X W (col) = 0.0d0
X Wi (col) = 0
X if (db) print 98, Fi (npivs+nfill),
X 1 Ipivcol (Fi (npivs+nfill)), nfill, Fx (npivs+nfill)
X endif
X 180 continue
X 170 continue
X
X do 109 k = 1, Size
X if (W (k) .ne. 0.0d0 .or. Wi (k) .ne. 0) then
X print *, 'ERROR in W', k, W (k), Wi (k)
X call crwait
X return
X endif
X 109 continue
X
Xc append fillin elements to R
X
X if (db) print *, 'storing fillin: nfill=', nfill
X do 200 k = 1, nfill
X lastp = sp
X call nextrow (sp, sk)
X if (sp .eq. Nil) then
X need = nfill - k
X needblocks = max (1, 1 + ((need-1) / Bl))
X if (db) then
X print *, 'allocating new space for fillin '
X print *, 'need', need, ' needblocks', needblocks
X endif
X memsize = needblocks * Rowbl
X MemTail = MemTail - memsize
X sp = MemTail
X if (MemTail .lt. MemHead) then
X print *, 'Out of row memory!', MemHead, MemTail
X return
X endif
X Row_next (lastp) = sp
X lastblock = lastp
X p = sp
X do 230 b = 1, needblocks
X nextblock = p + Rowbl
X Row_last (p) = lastblock
X Row_next (p) = nextblock
X lastblock = p
X 230 p = nextblock
X Row_next (lastblock) = Nil
X endif
X
X j = npivs + k
X Row_index (sp, sk) = Fi (j)
X Row_value (sp, sk) = Fx (j)
X if (db) then
X print *,k,sp,sk,'(',Fi(j),Ipivcol (Fi (j)), ')', Fx (j)
X endif
X 200 continue
X
X RowLen (row) = len + nfill
X if (db) call dumprow (row)
X 100 continue
X
Xc-------------------------------------------------------------------------------
Xc for each col to modify: (parallel loop)
X
X do 400 i = 1, NRedcol
X
X col = Redcol (i)
X p = ColHead (col)
X len = ColLen (col)
X
Xc scatter pivot columns into Wi
X
X if (db) then
X print *,'-----------------------------------------------'
X print *,'Reducing col ',col,' at i=',i,' p=',p,' len=',len
X print *,'Ipivcol (col)=', Ipivcol (col)
X print *, 'scatter C into Wi:'
X endif
X
X pivlink = Nil
X npivs = 0
X kl = 0
X do 410 j = 1, len, Bl
X kl = min (Bl, len - j + 1)
X do 420 k = 1, kl
X row = Col_index (p, k)
X piv = Ipivrow (row)
X if (piv .le. Plast) then
X Col_index (p, k) = pivlink
X pivlink = Address_of_Col_index (p, k)
X pl = Lp (piv)
X plen = Llen (piv)
X npivs = npivs + 1
X Fi (npivs) = piv
X if (db) then
X print *,'adding ',npivs,' pivot col ',row, ' ',piv
X print*,'plen=',plen,' pl=',pl,' Fi(npivs)=',Fi(npivs)
X endif
X do 430 l = 1, plen
X if (db) print 497, j, pl, l,
X 1 L_index (pl,l), Ipivrow (L_index (pl,l))
X 430 Wi (L_index (pl, l)) = 1
X if (db) print *, 'done adding pivot column ',row,piv
X Wi (row) = 0
X endif
X 420 continue
X lastp = p
X 410 p = Col_next (p)
X
X if (db) print *, 'number of pivot columns added to C:', npivs
X
Xc eliminate elements in pivot rows from C
X
X p = lastp
X k = kl
X 415 nextlink = Cm (pivlink)
X Cm (pivlink) = Col_index (p, k)
X pivlink = nextlink
X if (pivlink .eq. Nil) goto 425
X k = k - 1
X if (k .le. 0) then
X k = Bl
X p = idint (Col_last (p))
X endif
X goto 415
X 425 len = len - npivs
X sp = p
X sk = k
X ColLen (col) = len
X
Xc zero those nonzeros that are in C and Wi
X
X p = ColHead (col)
X do 450 j = 1, len, Bl
X kl = min (Bl, len - j + 1)
X do 460 k = 1, kl
X 460 Wi (Col_index (p, k)) = 0
X 450 p = Col_next (p)
X
Xc gather nonzeros still in Wi into Fi
X
X nfill = 0
X
X do 470 npiv = 1, npivs
X piv = Fi (npiv)
X pl = Lp (piv)
X plen = Llen (piv)
X if (db) then
X print *, ' '
X print *, 'fillin piv=',piv,' pl=',pl,' plen=',plen
X endif
X
X do 480 k = 1, plen
X row = L_index (pl, k)
X if (Wi (row) .ne. 0) then
X nfill = nfill + 1
X Fi (npivs + nfill) = row
X Wi (row) = 0
X if (db) print 498, Fi (npivs+nfill),
X 1 Ipivcol (Fi (npivs+nfill)), nfill
X endif
X 480 continue
X if (db) print *, 'done with fillin', piv
X 470 continue
X
X do 409 k = 1, Size
X if (Wi (k) .ne. 0) then
X print *, 'ERROR in Wi', k, Wi (k)
X call crwait
X return
X endif
X 409 continue
X
Xc append fillin elements to C
X
X if (db) print *, 'storing fillin: nfill=', nfill
X do 500 k = 1, nfill
X j = npivs + k
X Col_index (sp, sk) = Fi (j)
X
X if (db) print *,k,sp,sk,'(',Fi(j),Ipivrow(Fi(j)),')'
X lastp = sp
X call nextcol (sp, sk)
X
X if (sp .eq. Nil .and. k .ne. nfill) then
X need = nfill - k
X needblocks = max (1, 1 + ((need-1) / Bl))
X if (db) then
X print *, 'allocating new space for fillin '
X print *, 'need', need, ' needblocks', needblocks
X endif
X memsize = needblocks * Colbl
X MemCTail = MemCTail - memsize
X sp = MemCTail
X if (MemCTail .lt. MemCHead) then
X print *, 'Out of col memory!', MemCHead, MemCTail
X return
X endif
X Col_next (lastp) = sp
X lastblock = lastp
X p = sp
X do 530 b = 1, needblocks
X nextblock = p + Colbl
X Col_last (p) = lastblock
X Col_next (p) = nextblock
X lastblock = p
X 530 p = nextblock
X Col_next (lastblock) = Nil
X endif
X 500 continue
X
X ColLen (col) = len + nfill
X if (db) call dumpcol (col)
X 400 continue
X
Xc-------------------------------------------------------------------------------
Xc update Markowitz data structures
X
X do 600 i = 1, NRedrow
X row = Redrow (i)
X MRowflag (row) = 0
X nz = RowLen (row)
X p = MRowhead (nz)
X if (p .eq. Nil) then
X MRownext (row) = Nil
X MRowlast (row) = Nil
X MRowhead (nz) = row
X else
X MRownext (row) = p
X MRowlast (p) = row
X MRowlast (row) = Nil
X MRowhead (nz) = row
X endif
X if (nz .lt. MRowmin .or. MRowmin .eq. 0) MRowmin = nz
X 600 continue
X
X do 700 i = 1, NRedcol
X col = Redcol (i)
X MColflag (col) = 0
X nz = ColLen (col)
X p = MColhead (nz)
X if (p .eq. Nil) then
X MColnext (col) = Nil
X MCollast (col) = Nil
X MColhead (nz) = col
X else
X MColnext (col) = p
X MCollast (p) = col
X MCollast (col) = Nil
X MColhead (nz) = col
X endif
X if (nz .lt. MColmin .or. MColmin .eq. 0) MColmin = nz
X 700 continue
X
X if (db) then
X print *, 'Updated Markowitz data structure:'
X call dumpnzlist
X call checkrowcol
X endif
X
Xc-------------------------------------------------------------------------------
X
X return
X 95 format(' ',i4,i7,i4,' : ',i4,' ',i4,' ',g25.14,
X 1 ' <- pivot ', i7, i5)
X 97 format(' ',i4,i7,i4,' : ',i4,' ',i4,' ',g25.14)
X 99 format(' ',i4,1x,i4, ' ',g14.4,
X 1 ' = ',g14.4,' ',i1,' - ',g14.4,' * ',g14.4)
X 98 format(' ',i4,' ',i4,':',i4,' ',g25.14)
X
X 497 format(' ',i4,i7,i4,' : ',i4,' ',i4)
X 496 format(' ',i4,' (',i7,i4,') (',i7,i4,') ',i4,' ',i4)
X 498 format(' ',i4,' ',i4,':',i4)
X end
SHAR_EOF
true || echo 'restore of m4.F failed'
fi
# ============= main.F ==============
if test -f 'main.F' -a X"$1" != X"-c"; then
echo 'x - skipping main.F (File already exists)'
else
echo 'x - extracting main.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'main.F' &&
X call cmain
X end
SHAR_EOF
true || echo 'restore of main.F failed'
fi
# ============= marko.F ==============
if test -f 'marko.F' -a X"$1" != X"-c"; then
echo 'x - skipping marko.F (File already exists)'
else
echo 'x - extracting marko.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'marko.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
Xc markowitz: find the best pivot according to the markowitz criterion, and
Xc return the cost, and the pivot row and column
X
X subroutine markowitz (cost, prow, pcol)
X
X include 'global.h'
X integer cost, prow, pcol, row, col, kcost, i, nz, j
X integer p, minnz, pp, kk, k, ik, ip
X real*8 amax, au, findrowmax
X logical pr
X
X pr = Debug .eq. 1
X cost = Size ** 2
X p = Stage + NPivots
X minnz = min (MRowmin, MColmin)
X if (minnz .eq. 0) then
X cost = -1
X return
X endif
X
X if (pr) then
X print *, 'markowitz: Stage', Stage, ' NPivots', NPivots
X print *, ' Stage+NPivots:', p
X print *, ' minnz:', minnz
X print *, ' MColmin:',MColmin,' MRowmin:',MRowmin
X endif
X
X do 330 nz = minnz, Size
X
X if (pr) print *, 'searching nz=', nz
X
X if (cost .le. (nz - 1) ** 2) then
X if (pr) print *, 'returning after column search', cost
X return
X endif
X
Xc scan rows with nz nonzeros
X
X if (pr) print *, 'row scan, nz=', nz
X row = MRowhead (nz)
X10 if (row .eq. Nil) goto 30
X
Xc scan row for possible pivots
X
X amax = MRowmax (row)
X au = amax * 0.01d0
X k = 1
X p = RowHead (row)
X
X do 260 i = 1, RowLen (row)
X col = Row_index (p, k)
X kcost = (nz - 1) * (ColLen (col) - 1)
X
X if (MColflag (col) .eq. 1) then
X if (pr) print *, 'incompatible row pivot:', row, col
X
X elseif (kcost .lt. cost) then
X
X if (amax .lt. 0.0d0) then
X amax = findrowmax (row, col, ik, ip)
X MRowmax (row) = amax
X au = amax * 0.01d0
X endif
X
X if (abs (Row_value (p, k)) .gt. au) then
X cost = kcost
X prow = row
X pcol = col
X if (pr) print *, 'good row pivot:', prow,pcol,cost
X if (cost .le. (nz - 1) ** 2) then
X if (pr) then
X print *, 'final row pivot found'
X print *, 'amax:', amax, ' au:', au
X endif
X return
X endif
X elseif (pr) then
X print *, 'row pivot unstable:', row, col, kcost
X print *, Row_value (p, k), amax
X endif
X elseif (pr) then
X print *, 'row pivot too costly:', row, col, kcost
X endif
X call nextrow (p, k)
X260 continue
X
X25 row = MRownext (row)
X goto 10
X30 continue
X
X if (cost .le. nz * (nz - 1)) then
X if (pr) print *, 'returning after row search', cost
X return
X endif
X
Xc scan columns with nz nonzeros
X
X if (pr) print *, 'col scan, nz=', nz
X col = MColhead (nz)
X50 if (col .eq. Nil) goto 80
X
Xc scan column for possible pivots
X
X k = 1
X p = ColHead (col)
X
X do 310 i = 1, ColLen (col)
X row = Col_index (p, k)
X kcost = (nz - 1) * (RowLen (row) - 1)
X if (pr) print *, 'testing',row,col,' with cost ',kcost
X
X if (MRowflag (row) .eq. 1) then
X if (pr) print *, 'incompatible col pivot:', row, col
X
X elseif (kcost .lt. cost) then
X
X amax = MRowmax (row)
X if (amax .lt. 0.0d0) then
X amax = findrowmax (row, col, ik, ip)
X MRowmax (row) = amax
X else
X kk = 1
X pp = RowHead (row)
X do 305 j = 1, RowLen (row)
X if (col .eq. Row_index (pp, kk)) then
X ik = kk
X ip = pp
X goto 307
X endif
X call nextrow (pp, kk)
X305 continue
X endif
X307 au = amax * 0.01d0
X
X if (abs (Row_value (ip, ik)) .gt. au) then
X cost = kcost
X prow = row
X pcol = col
X if (pr) print *, 'good col pivot:', prow,pcol,cost
X if (cost .le. nz * (nz - 1)) then
X if (pr) then
X print *, 'final row pivot found'
X print *, 'amax:', amax, ' au:', au
X endif
X return
X endif
X elseif (pr) then
X print *, 'column pivot unstable:', row,col,kcost
X print *, Row_value (ip, ik), amax
X endif
X elseif (pr) then
X print *, 'column pivot too costly:', row, col, kcost
X endif
X call nextcol (p, k)
X310 continue
X
X col = MColnext (col)
X goto 50
X80 continue
X
X330 continue
X
X cost = -1
X return
X end
X
Xc-------------------------------------------------------------------------------
Xc pivotfound: remove the pivot row and column from the list of active row and
Xc columns (MRowhead and MColhead), and remove all rows with nonzeros in the
Xc pivot column and all columns with nonzeros in the pivot row.
X
X subroutine pivotfound (prow, pcol, nrow, ncol, whoami)
X include 'global.h'
X integer len, i, col, row, prow, pcol, nz, p, k, nrow, ncol
X integer whoami
Xc real*8 findrowmax, amax
X logical pr
X
X pr = Debug .eq. 1
X if (pr) print *, ' ======= Pivot found',prow,pcol,NPivots
X
X nrow = MRownext (prow)
X ncol = MColnext (pcol)
X
X k = 1
X p = RowHead (prow)
X do 20 i = 1, RowLen (prow)
X col = Row_index (p, k)
X if (MColflag (col) .eq. 0) then
X MColflag (col) = 1
X if (col .eq. ncol) ncol = MColnext (col)
X if (col .ne. pcol) then
X NRedcol = NRedcol + 1
X Redcol (NRedcol) = col
X endif
X len = ColLen (col)
X if (col .eq. MColhead (len)) then
X MColhead (len) = MColnext (col)
X MCollast (MColhead (len)) = Nil
X if (MColhead(len) .eq. Nil .and. len .eq. MColmin) then
X do 10 nz = MColmin+1, Size
X if (MColhead (nz) .ne. Nil) then
X MColmin = nz
X goto 15
X endif
X10 continue
X MColmin = 0
X endif
X else
X MColnext (MCollast (col)) = MColnext (col)
X MCollast (MColnext (col)) = MCollast (col)
X endif
X elseif (col .eq. pcol) then
X print *, 'ERROR: removing flagged pivot col'
X print *, prow, pcol, NPivots
X call crwait
X return
X endif
X15 continue
X call nextrow (p, k)
X20 continue
X
X k = 1
X p = ColHead (pcol)
X do 40 i = 1, ColLen (pcol)
X row = Col_index (p, k)
X if (MRowflag (row) .eq. 0) then
X if (row .eq. nrow) nrow = MRownext (row)
X if (row .ne. prow) then
X NRedrow = NRedrow + 1
X Redrow (NRedrow) = row
X endif
X MRowflag (row) = 1
X len = RowLen (row)
X if (row .eq. MRowhead (len)) then
X MRowhead (len) = MRownext (row)
X MRowlast (MRowhead (len)) = Nil
X if (MRowhead(len) .eq. Nil .and. len .eq. MRowmin) then
X do 30 nz = MRowmin+1, Size
X if (MRowhead (nz) .ne. Nil) then
X MRowmin = nz
X goto 35
X endif
X30 continue
X MRowmin = 0
X endif
X else
X MRownext (MRowlast (row)) = MRownext (row)
X MRowlast (MRownext (row)) = MRowlast (row)
X endif
X elseif (row .eq. prow) then
X print *, 'ERROR: removing flagged pivot row'
X print *, prow, pcol, NPivots
X call crwait
X return
X endif
X35 continue
X call nextcol (p, k)
X40 continue
X
X NPivots = NPivots + 1
X Prows (NPivots) = prow
X Pcols (NPivots) = pcol
X if (pr) then
X call dumpnzlist
X print *, ' ======= done with Pivot found',prow,pcol
X endif
X return
X end
SHAR_EOF
true || echo 'restore of marko.F failed'
fi
# ============= mattool.h ==============
if test -f 'mattool.h' -a X"$1" != X"-c"; then
echo 'x - skipping mattool.h (File already exists)'
else
echo 'x - extracting mattool.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'mattool.h' &&
X/*------------------------------------------------------------------------------
X SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
X This is freeware: you may use this tool and distribute it, as long as it's
X free. You can contact the author at na.tdavis@na-net.ornl.gov
X------------------------------------------------------------------------------*/
X#include <stdio.h>
X#include <strings.h>
X#include <pwd.h>
X#include <math.h>
X#include <ctype.h>
X#include <sys/types.h>
X#include <sys/dir.h>
X#include <sys/param.h>
X
X#include <suntool/sunview.h>
X#include <suntool/panel.h>
X#include <suntool/canvas.h>
X#include <suntool/tty.h>
X
X#define Filelen 512
X#define Line_len 200
X#define Maxfiles 200
X
X#define None 0
X#define Loaded 1
X#define Running 2
X#define Done 3
X
X#define StartDir "SparseMat"
X
X#define Until(condition) while (! (condition))
X
X#define Strncpy(dest,src,len) { strncpy (dest, src, len) ; dest [len] = '\0' ; }
X#define Malloc(size,type) \
X ((type *) malloc ((unsigned) (((size)+1) * (sizeof (type)))))
X#define Free(p) if (p) free ((char *) p)
X
X#define Button(label, proc, x, y) \
X panel_create_item (control_panel, PANEL_BUTTON, \
X PANEL_LABEL_IMAGE, panel_button_image (control_panel, label, 5,0), \
X PANEL_ITEM_X, ATTR_COL(x), PANEL_ITEM_Y, ATTR_ROW(y), \
X PANEL_NOTIFY_PROC, proc, 0)
X
X#define Slider(label, proc, value, min, max, y) \
X panel_create_item (control_panel, PANEL_SLIDER, \
X PANEL_LABEL_STRING, label, \
X PANEL_VALUE, value, PANEL_MIN_VALUE, min, PANEL_MAX_VALUE, max, \
X PANEL_ITEM_X, ATTR_COL(0), PANEL_ITEM_Y, ATTR_ROW(y), \
X PANEL_NOTIFY_PROC, proc, 0)
X
X#define Toggle(label, proc, value, col, row) \
X panel_create_item (control_panel, PANEL_TOGGLE, \
X PANEL_LABEL_STRING, label, PANEL_CHOICE_STRINGS, " ", 0, \
X PANEL_NOTIFY_PROC, proc, PANEL_TOGGLE_VALUE, 0, value, \
X PANEL_ITEM_X, ATTR_COL(col), PANEL_ITEM_Y, ATTR_ROW (row), 0) \
X
X#define Toggle_bit_on(value, bit) ((value) & (1 << (bit)))
X
X#define Text(label, proc, length, value, x, y) \
X panel_create_item (control_panel, PANEL_TEXT, \
X PANEL_VALUE_DISPLAY_LENGTH, length, \
X PANEL_LABEL_STRING, label, PANEL_VALUE, value, \
X PANEL_ITEM_X, ATTR_COL(x), PANEL_ITEM_Y, ATTR_ROW(y), \
X PANEL_NOTIFY_PROC, proc, 0)
X
X#define Set_msg(msg, string) panel_set (msg, PANEL_LABEL_STRING, string, 0)
X#define Clear_msg(msg) panel_set (msg, PANEL_LABEL_STRING, "\0", 0)
X#define Turn_off(item) panel_set (item, PANEL_SHOW_ITEM, FALSE, 0)
X#define Turn_on(item) panel_set (item, PANEL_SHOW_ITEM, TRUE, 0)
X
X#define Scrx(col) (2 + ((col) - 1) * ratio)
X#define Scry(row) (2 + ((row) - 1) * ratio)
X
X#define Set(x,y,xh,yh) pw_write (pw, x, y, xh, yh, PIX_SET, NULL, 0, 0)
X#define Clear(x,y,xh,yh) pw_write (pw, x, y, xh, yh, PIX_CLR, NULL, 0, 0)
X
X#define Set_element(row,col) Set (Scrx (col), Scry (row), elsize, elsize)
X#define Clear_element(row,col) Clear(Scrx (col), Scry (row), elsize, ratio)
X
X#define Set_row(row) Set (Scrx (1), Scry (row), Size*ratio, ratio)
X#define Clear_row(row) Clear(Scrx (1), Scry (row), Size*ratio, ratio)
X
X#define Gridrow(row) Set (Scrx(1), Scry(row)+ratio-1, Size*ratio, 1)
X#define Gridcol(col) Set (Scrx(col)+ratio-1, Scry(1), 1, Size*ratio)
X
X#define Set_rows(row,nrows) \
X Set (Scrx (1), Scry (row), Size*ratio, (nrows)*ratio)
X#define Clear_rows(row,nrows) \
X Clear (Scrx (1), Scry (row), Size*ratio, (nrows)*ratio)
X
X#define Clear_cols(col, ncols) \
X Clear (Scrx (col), Scry (1), (ncols)*ratio, Size*ratio)
X
X#define Compare(str1, str2) (strcmp (str1, str2) == 0)
X#define NCompare(str1, str2, n) (strncmp (str1, str2, n) == 0)
X#define Open_file(fp, name, mode) (fp = fopen (name, mode))
X#define Close_file(fp) (fclose (fp))
X#define Min(a,b) ((a) < (b) ? (a) : (b))
X#define Max(a,b) ((a) > (b) ? (a) : (b))
X#define Abs(a) ((a) > 0 ? (a) : -(a))
X#define Cost(i) ((RowLen(Pivrow(Stage+i))-1) * (ColLen(Pivcol(Stage+i))-1))
X
X/*----------------------------------------------------------------------------*/
X
Xdouble drand48 (), dget_value () ;
Xchar *malloc () ;
Xint quit_proc (), init_file (), fast_proc (), dir_proc (), file_proc (),
X dfile_proc (), null_proc (), get_value (), print_proc (),
X multi_proc (), mstep_proc (),
X clear_all (), display_multi_matrix (), display_all_multi_matrix ()
X;
X/*----------------------------------------------------------------------------*/
X
Xextern Frame frame ;
Xextern Panel control_panel ;
Xextern Canvas canvas ;
Xextern Panel_item dir_item, file_item, factor_item, dense_item ;
Xextern Menu file_menu ;
Xextern int wr, wc, wsize, n, ratio, elsize ;
Xextern Pixwin *pw ;
Xextern FILE *matfile ;
Xextern char mat_title [80], key [10], type [10], path [256] ;
Xextern int i, col1, col2, row1, row2,
X cols, nonzeros, rhs, fail, drops, printing, ops, init_nz, print_level ;
X
Xextern DIR *cwd ;
Xextern char pathname [MAXPATHLEN], filename [32] ;
Xextern char files [Maxfiles][32], smsize [Maxfiles][8] ;
Xextern int entries, msize [Maxfiles], status, display ;
Xextern int mstep ;
X
X/*----------------------------------------------------------------------------*/
X
X#include "cdefine.h"
SHAR_EOF
true || echo 'restore of mattool.h failed'
fi
# ============= mback.F ==============
if test -f 'mback.F' -a X"$1" != X"-c"; then
echo 'x - skipping mback.F (File already exists)'
else
echo 'x - extracting mback.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'mback.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
X subroutine mdone
X include 'global.h'
X call mback (Xm (Densep), Densesize, Sparsesize)
X return
X end
X
X subroutine mback (a, n, s)
X include 'global.h'
X integer row, col, p, len, k, n, s
Xc integer i, kl, pl, pu, ku
Xc real*8 maxz
X real*8 y (Nmax), xtmax, a (n, n)
X logical pr
X equivalence (y, W)
X pr = Printflag .eq. 1
X
Xc LU factorization complete, reorder indices of sparse L and U
X
X do 10 col = 1, Sparsesize
X p = Lp (col)
X len = Llen (col)
X do 10 k = 1, len
X 10 L_index (p, k) = Ipivrow (L_index (p, k))
X
X do 20 row = 1, Sparsesize
X p = Up (row)
X len = Ulen (row)
X do 20 k = 1, len
X 20 U_index (p, k) = Ipivcol (U_index (p, k))
X
Xcc find Z = LU = PAQ
Xc
Xc if (pr .and. Size .le. 8) then
Xc do 21 row = 1, Size
Xc do 21 col = 1, Size
Xc 21 Z (row, col) = 0.0d0
Xc
Xc do 22 i = 1, Size
Xc pu = Up (i)
Xc pl = Lp (i)
Xc do 22 ku = 1, Ulen (i)
Xc col = U_index (pu, ku)
Xc row = i
Xc Z (row, col) = Z (row, col) + U_value (pu, ku)
Xc do 22 kl = 1, Llen (i)
Xc row = L_index (pl, kl)
Xc 22 Z(row,col)=Z(row,col)+L_value(pl,kl)*U_value(pu,ku)
Xc
Xc print *, 'L times U:'
Xc call dumpz
Xc
Xc do 24 k = 1, Nonzeros
Xc row = Ipivrow (A_row (k))
Xc col = Ipivcol (A_col (k))
Xc 24 Z (row, col) = Z (row, col) - A_value (k)
Xc
Xc print *, 'PAQ - LU: (should be close to zero)'
Xc call dumpz
Xc
Xc maxz = 0.0d0
Xc do 25 row = 1, Size
Xc do 25 col = 1, Size
Xc maxz = max (maxz, abs (Z (row, col)))
Xc 25 Z (row, col) = 0.0d0
Xc print *, 'largest value in (PAQ-LU) = ', maxz
Xc
Xc do 26 k = 1, Nonzeros
Xc row = Ipivrow (A_row (k))
Xc col = Ipivcol (A_col (k))
Xc 26 Z (row, col) = A_value (k)
Xc
Xc print *, 'PAQ:'
Xc call dumpz
Xc call dumpperm
Xc endif
X
Xc solve the lower triangular system (Ly = Pb) for y
X
Xc print *, 'dumpdense in mback: s=',s, ' n=',n
Xc call dumpdense (a, n, s, Pivrow, Ipivrow)
X
X if (pr) print *, 'Pb:'
X do 30 row = 1, Size
X y (row) = Rhs (Pivrow (row))
X 30 if (pr) print *, row, Pivrow (row), y (row)
X
X if (pr) print *, 'solve Ly = Pb for y:'
X
X do 40 col = 1, Sparsesize
X p = Lp (col)
X len = Llen (col)
X if (pr) print *, col, y (col)
X do 40 k = 1, len
X row = L_index (p, k)
X 40 y (row) = y (row) - y (col) * L_value (p, k)
X
X do 45 col = Sparsesize + 1, Size
X if (pr) print *, col, y (col)
X do 45 row = col+1, Size
X 45 y (row) = y (row) - y (col) * a (row-s, col-s)
X
Xc solve the upper triangular system (Uz = y) for z, then x = Qz
X
X if (pr) print *, 'solve Uz = y for z:', Size, Sparsesize
X
X do 55 row = Size, Sparsesize + 1, -1
X do 65 col = row+1, Size
X 65 y (row)= y (row) - y (col) * a (row-s, col-s)
X y (row) = y (row) / a (row-s, row-s)
X 55 if (pr) print *, row, y (row), a (row-s,row-s)
X
X do 60 row = Sparsesize, 1, -1
X p = Up (row)
X len = Ulen (row)
X do 50 k = 2, len
X col = U_index (p, k)
X 50 y (row) = y (row) - y (col) * U_value (p, k)
X y (row) = y (row) / U_value (p, 1)
X 60 if (pr) print *, row, y (row)
X
X if (pr) print *, 'x = Qz:'
X do 70 row = 1, Size
X Rhs (Pivcol (row)) = y (row)
X 70 if (pr) print *, row, Pivcol (row), y (row)
X
Xc check the accuracy of the solution
X
X Error = 0.0d0
X xtmax = 0.0d0
X do 80 row = 1, Size
X xtmax = max (xtmax, abs (Xt (row)))
X Error = max (Error, abs (Rhs (row) - Xt (row)))
X 80 if (pr) print *, Xt (row), Rhs (row), Error, xtmax
X Error = Error / xtmax
X return
X end
SHAR_EOF
true || echo 'restore of mback.F failed'
fi
# ============= read.F ==============
if test -f 'read.F' -a X"$1" != X"-c"; then
echo 'x - skipping read.F (File already exists)'
else
echo 'x - extracting read.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'read.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
Xc-----------------------------------------------------------------------
Xc read: reads a matrix from the Harwell/Boeing collection into im/xm.
Xc Does not assume any ordering, and does not give an ordered output.
Xc
Xc Format in cm/im/xm:
Xc rk cm (k) A_row (k)
Xc ck im (k) A_col (k)
Xc A (rk,ck) xm (k) A_value (k)
Xc
Xc-----------------------------------------------------------------------
X
X subroutine readmatrix ()
X
X include 'global.h'
X
X character title*72, type*3, ptrfmt*16, indfmt*16, valfmt*20,
X 1 rhsfmt*20
X integer totcrd, ptrcrd, indcrd, valcrd, rhscrd, nrow,
X 1 nrhs, i, col, k, colptr (Nmax), drops
X logical valued
X real*8 drand, skew, x
X equivalence (colptr, Iw1)
X
X print *, 'Reading matrix...'
X open (unit = 7, file = "matrix.in", status = "old")
X read (7, 40)
X 1 title, Name,
X 1 totcrd, ptrcrd, indcrd, valcrd, rhscrd,
X 1 type, nrow, Size, Nonzeros, nrhs,
X 1 ptrfmt, indfmt, valfmt, rhsfmt
X
X valued = valcrd .gt. 0
X if (nrow.ne.Size .or. Size.gt.Nmax .or. Nonzeros.gt.Mmax) then
X print *, 'Invalid matrix -- rectangular or too big'
X Size = - Size
X close (unit = 7)
X return
X endif
X
Xc allocate memory for original A matrix at top of free memory
X
X MemTail = Mmax+1
X MemHead = Nonzeros+1
X MemCTail = Cmax+1
X MemCHead = Nonzeros+1
X
Xc read column pointers, row indices, and nonzero values
X
X read (7, ptrfmt) (colptr (col), col = 1, Size+1)
X read (7, indfmt) (A_row (k), k = 1, Nonzeros)
X if (valued) then
X read (7, valfmt) (A_value (k), k = 1, Nonzeros)
X else
X print *, 'Creating random values for matrix pattern'
Xc restart random-number generator
X x = drand (1)
X do 10 k = 1, Nonzeros
X 10 A_value (k) = drand (0) + 0.01d0
X endif
X close (unit = 7)
X
Xc find column indices
X
Xcvd$l nodepchk
X do 20 col = 1, Size
X do 20 i = colptr (col), colptr (col+1) - 1
X 20 A_col (i) = col
X
Xc delete zero elements
X
X drops = 0
X do 30 k = 1, Nonzeros
X if (A_value (k) .eq. 0.0d0) then
X drops = drops + 1
X elseif (drops .gt. 0) then
X A_row ((k-drops)) = A_row (k)
X A_col ((k-drops)) = A_col (k)
X A_value ((k-drops)) = A_value (k)
X endif
X 30 continue
X Nonzeros = Nonzeros - drops
X
Xc add upper triangular part if symmetric or skew symmetric
X
X skew = 0.0
X if (type (2:2) .eq. 'z' .or. type (2:2) .eq. 'Z') skew = -1.0
X if (type (2:2) .eq. 's' .or. type (2:2) .eq. 'S') skew = 1.0
X
X if (skew .ne. 0.0) then
X i = Nonzeros
X do 60 k = 1, Nonzeros
X if (A_row (k) .ne. A_col (k)) then
X i = i + 1
X A_row (i) = A_col (k)
X A_col (i) = A_row (k)
X A_value (i) = skew * A_value (k)
X endif
X60 continue
X Nonzeros = i
X MemHead = Nonzeros + 1
X MemCHead = Nonzeros + 1
X endif
X
Xc print results
X
X print 50, title (2:72), Name, type, Size, Nonzeros, valued
X return
X
X 40 format (a72, a8 / 5i14 / a3, 11x, 4i14 / 2a16, 2a20)
X 50 format (' ', a71, a8 / ' type:', a3, ' size:', i6,
X 1 ' nonzeros:', i8, ' valued:', i2)
X end
SHAR_EOF
true || echo 'restore of read.F failed'
fi
# ============= store.F ==============
if test -f 'store.F' -a X"$1" != X"-c"; then
echo 'x - skipping store.F (File already exists)'
else
echo 'x - extracting store.F (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'store.F' &&
XC-------------------------------------------------------------------------------
XC SparseDyanmics tool, copyright 1990 by Tim Davis, all rights reserved.
XC This is freeware: you may use this tool and distribute it, as long as it's
XC free. You can contact the author at na.tdavis@na-net.ornl.gov
XC-------------------------------------------------------------------------------
X#include "fdefine.h"
Xc-------------------------------------------------------------------------------
Xc storematrix: store the A matrix from the front of the memory space to
Xc the link list data structure. The A matrix is not in any particular
Xc order
X
X subroutine storematrix
X
X include "global.h"
X integer k, i, p, b, nblocks, lastblock, nextblock, row, col, nz
X integer ci (Nmax+2), cp (Nmax+2), memsize
X integer ri (Nmax+2), rp (Nmax+2)
X real*8 x, value
X equivalence (ci, Iw1), (cp, Iw2), (ri, Iw3), (rp, Iw4)
X
Xc initialize data structures
X
X x = 1.0 / dble (2 * Size)
X NStages = 0
X Stage = 0
X Densesize = 0
X Sparsesize = Size
X Densep = 0
X
X if (Size .le. 8) then
X do 5 i = 1, Size
X 5 Xt (i) = i
X else
X do 6 i = 1, Size
X 6 Xt (i) = 1.0d0 + (i * 0.01d0)
X endif
X
X do 10 i = 1, Size
X Wi (i) = 0
X W (i) = 0.0
X Rhs (i) = 0.0d0
X RowLen (i) = 0
X ColLen (i) = 0
X Pivrow (i) = i
X Pivcol (i) = i
X Ipivrow (i) = i
X Ipivcol (i) = i
X MRowhead (i) = Nil
X MRownext (i) = Nil
X MRowlast (i) = Nil
X MRowflag (i) = 0
X MRowmax (i) = -1.0d0
X MColhead (i) = Nil
X MColnext (i) = Nil
X MCollast (i) = Nil
X MColflag (i) = 0
X ci (i) = 0
X 10 ri (i) = 0
X
Xc find the number of nonzeros in each row and column of A, find Rhs
X
X do 20 k = 1, Nonzeros
X Rhs (A_row (k)) = Rhs (A_row (k))
X 1 + (Xt (A_col (k)) * A_value (k))
X RowLen (A_row (k)) = RowLen (A_row (k)) + 1
X 20 ColLen (A_col (k)) = ColLen (A_col (k)) + 1
X if (Debug .eq. 1) call dumprhs
X
Xc initialize the Markowitz selection data structure
X
X do 22 row = 1, Size
X nz = RowLen (row)
X p = MRowhead (nz)
X if (p .eq. Nil) then
X MRownext (row) = Nil
X MRowlast (row) = Nil
X MRowhead (nz) = row
X else
X MRownext (row) = p
X MRowlast (p) = row
X MRowlast (row) = Nil
X MRowhead (nz) = row
X endif
X 22 continue
X
X do 24 col = 1, Size
X nz = ColLen (col)
X p = MColhead (nz)
X if (p .eq. Nil) then
X MColnext (col) = Nil
X MCollast (col) = Nil
X MColhead (nz) = col
X else
X MColnext (col) = p
X MCollast (p) = col
X MCollast (col) = Nil
X MColhead (nz) = col
X endif
X 24 continue
X
X do 26 nz = 1, Size
X if (MRowhead (nz) .ne. 0) then
X MRowmin = nz
X goto 27
X endif
X 26 continue
X MRowmin = 0
X
X 27 do 28 nz = 1, Size
X if (MColhead (nz) .ne. 0) then
X MColmin = nz
X goto 29
X endif
X 28 continue
X MColmin = 0
X 29 continue
X
Xc allocate memory blocks for the rows of A
X
X do 40 i = 1, Size
X nblocks = max (1, 1 + ((RowLen(i)-1) / Bl))
X memsize = nblocks * Rowbl
X MemTail = MemTail - memsize
X p = MemTail
X if (MemTail .lt. MemHead) then
X print *, 'Store A: Out of row memory!', MemHead, MemTail
X call crwait
X return
X endif
X RowHead (i) = p
X rp (i) = p
X lastblock = -i
X do 30 b = 1, nblocks
X nextblock = p + Rowbl
X Row_last (p) = lastblock
X Row_next (p) = nextblock
X lastblock = p
X 30 p = nextblock
X 40 Row_next (lastblock) = Nil
X
Xc allocate memory blocks for the columns of A
X
X do 60 i = 1, Size
X nblocks = max (1, 1 + ((ColLen(i)-1) / Bl))
X memsize = nblocks * Colbl
X MemCTail = MemCTail - memsize
X p = MemCTail
X if (MemCTail .lt. MemCHead) then
X print *, 'Store A: Out of col memory!', MemCHead, MemCTail
X call crwait
X return
X endif
X ColHead (i) = p
X cp (i) = p
X lastblock = -i
X do 50 b = 1, nblocks
X nextblock = p + Colbl
X Col_last (p) = lastblock
X Col_next (p) = nextblock
X lastblock = p
X 50 p = nextblock
X 60 Col_next (lastblock) = Nil
X
Xc store the nonzero elements from original structure to link list structure
X
X do 70 k = 1, Nonzeros
X row = A_row (k)
X col = A_col (k)
X value = A_value (k)
X
X ci (col) = ci (col) + 1
X if (ci (col) .gt. Bl) then
X ci (col) = 1
X cp (col) = Col_next (cp (col))
X if (cp (col) .lt. MemCTail) print *, 'cp (col) ERROR ****',
X 1 row, col, cp (col), ci (col)
X endif
X Col_index (cp (col), ci (col)) = row
X
X ri (row) = ri (row) + 1
X if (ri (row) .gt. Bl) then
X ri (row) = 1
X rp (row) = Row_next (rp (row))
X if (rp (row) .lt. MemTail) print *, 'rp (row) ERROR ****',
X 1 row, col, rp (row), ri (row)
X endif
X Row_index (rp (row), ri (row)) = col
X 70 Row_value (rp (row), ri (row)) = value
X
Xc free memory space used for the original matrix A
X
Xc MemHead = 1
Xc MemCHead = 1
X
X if (Printflag .eq. 1) then
X call dumpmatrix
X call dumpnzlist
X endif
X return
X end
SHAR_EOF
true || echo 'restore of store.F failed'
fi
exit 0