Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

Add function to write RB data to disk.

- Ability to write matrices and supplementary data in RB format.
- Now depends upon fortranformat (from PyPi) instead of
  FortranFormat (from ScientificPython).
  • Loading branch information...
commit c4120158ca5a4a4676639256d84f828f03ee338b 1 parent 17a6b80
Dominique authored

Showing 1 changed file with 312 additions and 34 deletions. Show diff stats Hide diff stats

  1. +312 34 pyorder/tools/hrb.py
346 pyorder/tools/hrb.py
@@ -45,6 +45,7 @@
45 45 """
46 46 import numpy
47 47 import FortranFormat as F
  48 +from fortranformat import FortranRecordReader, FortranRecordWriter
48 49
49 50 class HarwellBoeingMatrix:
50 51 """
@@ -83,8 +84,12 @@ def __init__(self, fname, **kwargs):
83 84 self.readMatrix(fp, **kwargs)
84 85 fp.close()
85 86
  87 + def data(self):
  88 + #if self.ip is None or self.ind is None: return (None,None,None)
  89 + return (self.ip, self.ind, self.val)
  90 +
86 91 def find(self):
87   - if self.ip is None or self.ind is None: return None
  92 + if self.ip is None or self.ind is None: return (None,None)
88 93 row = numpy.empty(self.nnzero, numpy.int)
89 94 col = numpy.empty(self.nnzero, numpy.int)
90 95 # Adjust indices to 0-based scheme
@@ -92,9 +97,9 @@ def find(self):
92 97 for j in range(self.ip[i], self.ip[i+1]):
93 98 col[j] = i
94 99 row[j] = self.ind[j]
95   - return(row,col)
  100 + return (row,col)
96 101
97   - def readArray(self, fp, which, nelm, format):
  102 + def _Old_readArray(self, fp, which, nelm, format):
98 103 fmt = F.FortranFormat(str(format))
99 104 ind = 0
100 105 while ind < nelm-1:
@@ -104,11 +109,31 @@ def readArray(self, fp, which, nelm, format):
104 109 ind = ind2
105 110 return
106 111
107   - def fortranRead(self, stream, format):
  112 + def _Old_fortranRead(self, stream, format):
108 113 fmt = F.FortranFormat(format)
109 114 fdata = F.FortranLine(stream, fmt)
110 115 return fdata.data
111 116
  117 + def readArray(self, fp, which, nelm, format):
  118 + #print 'Reading %d values with format %s' % (nelm, format)
  119 + fmt = FortranRecordReader(format)
  120 + ind = 0
  121 + while ind < nelm-1:
  122 + fdata = fmt.read(fp.readline())
  123 + ind2 = min(ind + len(fdata), nelm)
  124 + which[ind:ind2] = fdata[:ind2-ind]
  125 + ind = ind2
  126 + # Read last line, if any.
  127 + if ind < nelm:
  128 + fdata = fmt.read(fp.readline())
  129 + which[ind:] = fdata[:nelm-ind]
  130 + return
  131 +
  132 + def fortranRead(self, stream, format):
  133 + fmt = FortranRecordReader(format)
  134 + fdata = fmt.read(stream)
  135 + return fdata
  136 +
112 137 def readMatrix(self, fp, **kwargs):
113 138 self.patternOnly = kwargs.get('patternOnly', False)
114 139 self.readRhs = kwargs.get('readRhs', False)
@@ -267,23 +292,26 @@ def readMatrix(self, fp, **kwargs):
267 292 else:
268 293
269 294 # Read supplementary data
270   - self.dattyp, positn, orgniz, caseid, numerf, m, nvec, nauxd = \
271   - fRead(buffer1, 'A3, 2A1, 1X, A8, 2X, A1, 3(2X, I13)')
  295 + (self.dattyp, self.positn, self.orgniz, self.caseid, numerf, m,
  296 + nvec, self.nauxd) = fRead(buffer1,
  297 + 'A3,2A1,1X,A8,2X,A1,3(2X,I13)')
272 298 auxfm1, auxfm2, auxfm3 = fRead(buffer2, '3A20')
  299 + self.nrow = m
  300 + self.nvec = self.ncol = nvec
273 301
274 302 # Read integer data
275   - if (self.dattyp=='rhs' and orgniz=='s') or \
  303 + if (self.dattyp=='rhs' and self.orgniz=='s') or \
276 304 self.dattyp in ['ipt','icv','ord']:
277 305 if self.dattyp=='ord':
278   - nauxd = m * nvec
  306 + self.nauxd = m * nvec
279 307 auxfm = auxfm1
280 308 else:
281 309 self.ip = numpy.empty(nvec+1, numpy.int)
282 310 self.readArray(fp, self.ip, nvec+1, auxfm1)
283 311 auxfm = auxfm2
284 312 self.ip -= 1 # Adjust to 0-based indexing
285   - self.ind = numpy.empty(nauxd, numpy.int)
286   - self.readArray(fp, self.ind, nauxd, auxfm)
  313 + self.ind = numpy.empty(self.nauxd, numpy.int)
  314 + self.readArray(fp, self.ind, self.nauxd, auxfm)
287 315 self.ind -= 1 # Adjust to 0-based indexing
288 316 if self.dattyp != 'rhs': return
289 317
@@ -295,17 +323,257 @@ def readMatrix(self, fp, **kwargs):
295 323 dataType = numpy.complex
296 324 elif numerf=='i':
297 325 dataType = numpy.int
298   - if self.dattyp != 'rhs': nauxd = m*nvec
299   - if self.dattyp == 'rhs' and orgniz == 's':
  326 + if self.dattyp != 'rhs': self.nauxd = m*nvec
  327 + if self.dattyp == 'rhs' and self.orgniz == 's':
300 328 auxfm = auxfm3
301 329 else:
302 330 auxfm = auxfm1
303   - self.val = numpy.empty(nauxd, dataType)
304   - self.readArray(fp, self.val, nauxd, auxfm)
  331 + self.val = numpy.empty(self.nauxd, dataType)
  332 + self.readArray(fp, self.val, self.nauxd, auxfm)
  333 + self.nnzero = self.nauxd
305 334
306 335 return
307 336
308 337
  338 +# Helper functions.
  339 +
  340 +def get_int_fmt(n):
  341 + "Return appropriate format for integer arrays."
  342 + fmts = ['(40I2)', '(26I3)', '(20I4)', '(16I5)', '(13I6)', '(11I7)',
  343 + '(10I8)', '(8I9)', '(8I10)', '(7I11)', '(4I20)']
  344 + nlines = [40,26,20,16,13,11,10,8,8,7,4]
  345 + nn = n
  346 + for k in range(1,n+1):
  347 + if nn < 10: break
  348 + nn /= 10
  349 + if k <= 10: return (fmts[k+1], nlines[k+1])
  350 + return (fmts[10], nlines[10])
  351 +
  352 +def get_real_fmt(p):
  353 + "Return appropriate format for real array (1 <= p <= 17)."
  354 + fmt = ['(8E10.1E3)', '(7E11.2E3)', '(6E12.3E3)', '(6E13.4E3)',
  355 + '(5E14.5E3)', '(5E15.6E3)', '(5E16.7E3)', '(4E17.8E3)',
  356 + '(4E18.9E3)', '(4E19.10E3)', '(4E20.11E3)', '(3E21.12E3)',
  357 + '(3E22.13E3)', '(3E23.14E3)', '(3E24.15E3)', '(3E25.16E3)']
  358 + fmt1 = ['(1P,8E10.1E3)', '(1P,7E11.2E3)', '(1P,6E12.3E3)',
  359 + '(1P,6E13.4E3)', '(1P,5E14.5E3)', '(1P,5E15.6E3)',
  360 + '(1P,5E16.7E3)', '(1P,4E17.8E3)', '(1P,4E18.9E3)',
  361 + '(1P,4E19.10E3)', '(1P,4E20.11E3)', '(1P,3E21.12E3)',
  362 + '(1P,3E22.13E3)', '(1P,3E23.14E3)', '(1P,3E24.15E3)',
  363 + '(1P,3E25.16E3)']
  364 + lens = [8,7,6,6,5,5,5,4,4,4,4,3,3,3,3,3]
  365 + return (fmt[p-2], fmt1[p-2], lens[p-2])
  366 +
  367 +def fortranWriteLine(data, stream, fformat):
  368 + "Write `data` to `stream` according to Fortran format `fformat`."
  369 + fmt = FortranRecordWriter(fformat)
  370 + stream.write(fmt.write(data))
  371 + stream.write('\n')
  372 + return
  373 +
  374 +def fortranWriteArray(data, chunk_size, stream, fformat):
  375 + """
  376 + Write array `data` to `stream`, possibly using multiple lines,
  377 + according to Fortran format `fformat`.
  378 + """
  379 + nelts = len(data)
  380 + nelts_per_line = nelts/chunk_size
  381 + #print 'Writing %d elements %d per line...' % (nelts, chunk_size)
  382 + for k in range(nelts_per_line):
  383 + chunk = data[k*chunk_size:(k+1)*chunk_size]
  384 + fortranWriteLine(chunk, stream, fformat)
  385 + if nelts_per_line*chunk_size < nelts:
  386 + fortranWriteLine(data[nelts_per_line*chunk_size:], stream, fformat)
  387 + return
  388 +
  389 +# End of helper functions.
  390 +
  391 +
  392 +def write_rb(fname, nrow, ncol, ip, ind,
  393 + val=None, precision=17, symmetric=False, skew=False,
  394 + title='Generic', key='Generic'):
  395 + """
  396 + Write a sparse matrix to file in Rutherford-Boeing format. The input matrix
  397 + must be described in compressed column format by the arrays `ip` and `ind`.
  398 + Rows must be ordered in each column.
  399 + If numerical values are to be written to file, they should be specified in
  400 + the array `val`.
  401 +
  402 + Currently only supports assembled matrices in compressed column format.
  403 + """
  404 +
  405 + patternOnly = (val == None)
  406 + rectangular = (nrow != ncol)
  407 + ne = len(ind)
  408 +
  409 + # Check that columns are in order.
  410 +# for j in range(ncol):
  411 +# if ip[j+1] <= ip[j]:
  412 +# raise ValueError, 'Columns must be ordered.'
  413 +
  414 + # Check that rows are in order in each column.
  415 +# for j in range(ncol):
  416 +# for k in range(ip[j], ip[j+1]-1):
  417 +# if ind[k] >= ind[k+1]:
  418 +# raise ValueError, 'Rows must be ordered in each column.'
  419 +
  420 + # Set mxtype.
  421 + mxtype0 = 'r'
  422 + if patternOnly: mxtype0 = 'p'
  423 + if symmetric: mxtype1 = 's'
  424 + if skew: mxtype1 = 'z'
  425 + if rectangular: mxtype1 = 'r'
  426 +
  427 + mxtype = mxtype0 + mxtype1 + 'a'
  428 +
  429 + # Set format and number card images for pointer array.
  430 + (ptrfmt, ptrn) = get_int_fmt(ne+1)
  431 + ptrcrd = ncol/ptrn + 1
  432 +
  433 + # Set format and number card images for index array.
  434 + (indfmt, indn) = get_int_fmt(nrow)
  435 + indcrd = (ne-1)/indn + 1
  436 +
  437 + # Set number of card images for numerical entries.
  438 + if patternOnly:
  439 + valcrd = 0 ; valfmi = ' '
  440 + else:
  441 + (valfmi, valfmo, valn) = get_real_fmt(precision)
  442 + valcrd = (ne-1)/valn + 1
  443 +
  444 + totcrd = ptrcrd + indcrd + valcrd
  445 + neltvl = 0
  446 +
  447 + fp = open(fname, 'w')
  448 +
  449 + lt = len(title)
  450 + if lt < 72: title = title + (72-lt)*' '
  451 + lk = len(key)
  452 + if lk < 8: key = key + (8-lk)*' '
  453 +
  454 + # Write header.
  455 + fortranWriteLine([title, key], fp, 'A72,A8')
  456 + fortranWriteLine([totcrd, ptrcrd, indcrd, valcrd], fp, 'I14,3(1X,I13)')
  457 + fortranWriteLine([mxtype, nrow, ncol, ne, neltvl], fp, 'A3,11X,4(1X,I13)')
  458 + fortranWriteLine([ptrfmt, indfmt, valfmi], fp, '2A16,A20')
  459 +
  460 + # Write pointer and index arrays. Ajust for 1-based indexing.
  461 + #print 'Writing pointer array...'
  462 + fortranWriteArray(ip+1, ptrn, fp, ptrfmt)
  463 + #print 'Writing index array...'
  464 + fortranWriteArray(ind+1, indn, fp, indfmt)
  465 +
  466 + # Write matrix entries.
  467 + neltvl = ne
  468 + if not patternOnly:
  469 + #print 'Writing matrix entries...'
  470 + fortranWriteArray(val, valn, fp, valfmo)
  471 +
  472 + fp.close()
  473 + return
  474 +
  475 +
  476 +def write_aux(fname, nrow, nvec, precision=17, title='Generic', key='Generic',
  477 + caseid='Generic', dattyp='rhs', positn='r', orgniz='d', nauxd=None,
  478 + ip=None, ind=None, val=None):
  479 + """
  480 + Write supplementary data to file in Rutherford-Boeing format.
  481 +
  482 + Only real data is supported for now.
  483 + """
  484 +
  485 + data_types = ['ord','rhs','sln','est','evl','svl','evc','svc','sbv','sbm',
  486 + 'sbp','ipt','icv','lvl','geo','avl']
  487 + organizations = ['s','d','e']
  488 + positions = ['r','l','s']
  489 +
  490 + if dattyp not in data_types:
  491 + raise ValueError, 'Unknown data type: %s' % dattyp
  492 + if positn not in positions:
  493 + raise ValueError, 'Unknown position: %s' % positn
  494 + if orgniz not in organizations:
  495 + raise ValueError, 'Unknown organization: %s' % orgniz
  496 +
  497 + if dattyp in ['evl','svl','lvl','sbp']: nvec = 1
  498 + if dattyp in ['evl','svl','lvl','sbv','sbm','sbp','avl']: positn = ' '
  499 + if dattyp != 'rhs': orgniz = ' '
  500 +
  501 + numerf = 'r'
  502 + if dattyp == 'ord': numerf = 'i'
  503 + if dattyp in ['ipt','icv']: numerf = 'p'
  504 +
  505 + if orgniz != 'e':
  506 + nauxd = 0
  507 + if orgniz == 's' or dattyp in ['icv','ipt']:
  508 + if ip is None:
  509 + raise ValueError, 'Need pointer array for data type %s' % dattyp
  510 + nauxd = ip(nvec)-1
  511 + if orgniz == 'd': nauxd = nrow*nvec
  512 +
  513 + # Set data formats.
  514 + auxfm1 = auxfm2 = auxfm3 = ' '
  515 + fm1 = fm3 = ' '
  516 +
  517 + if dattyp in ['ipt','icv']:
  518 + (auxfm1, n1) = get_int_fmt(nauxd+1)
  519 + (auxfm2, n2) = get_int_fmt(nrow)
  520 + elif dattyp == 'ord':
  521 + (auxfm1, n1) = get_int_fmt(nrow)
  522 + else:
  523 + if precision < 2 or precision > 17: precision = 17
  524 + if dattyp == 'rhs' and orgniz == 's':
  525 + (auxfm1, n1) = get_int_fmt(nauxd+1)
  526 + (auxfm2, n2) = get_int_fmt(nrow)
  527 + (auxfm3, fm3, n3) = get_real_fmt(precision)
  528 + else:
  529 + (auxfm1, fm1, n1) = get_real_fmt(precision)
  530 +
  531 + fp = open(fname, 'w')
  532 +
  533 + lt = len(title)
  534 + if lt < 72: title = title + (72-lt)*' '
  535 + lk = len(key)
  536 + if lk < 8: key = key + (8-lk)*' '
  537 +
  538 + # Write header.
  539 + fortranWriteLine([title, key], fp, 'A72,A8')
  540 + fortranWriteLine([dattyp, positn, orgniz, caseid, numerf, nrow, nvec,
  541 + nauxd], fp, 'A3,2A1,1X,A8,1X,A1,3(1X,I13)')
  542 + fortranWriteLine([auxfm1, auxfm2, auxfm3], fp, '3A20')
  543 +
  544 + # Write pointer and index arrays. Ajust for 1-based indexing.
  545 + if (dattyp == 'rhs' and orgniz == 's') or dattyp in ['ipt','icv']:
  546 + #print 'Writing pointer array...'
  547 + fortranWriteArray(ip+1, n1, fp, auxfm1)
  548 + #print 'Writing index array...'
  549 + fortranWriteArray(ind+1, n2, fp, auxfm2)
  550 +
  551 + # Write entries.
  552 + #print 'Writing entries...'
  553 + if dattyp == 'rhs' and orgniz == 's':
  554 + fortranWriteArray(val, n3, fp, fm3)
  555 + else:
  556 + fortranWriteArray(val, n1, fp, fm1)
  557 +
  558 + fp.close()
  559 + return
  560 +
  561 +
  562 +def write_aux_from_rb(fname, rbdata):
  563 + (ip, ind, val) = rbdata.data()
  564 + write_aux(fname, rbdata.nrow, ip.shape[0], dattyp=rbdata.dattyp,
  565 + title=rbdata.title, key=rbdata.key, caseid=rbdata.caseid,
  566 + positn=rbdata.positn, orgniz=rbdata.orgniz, nauxd=rbdata.nauxd,
  567 + ip=ip, ind=ind, val=val)
  568 +
  569 +def write_aux_from_numpy(fname, array, **kwargs):
  570 + if len(array.shape) == 1:
  571 + nvec = 1
  572 + else:
  573 + nvec = array.shape[1]
  574 + write_aux(fname, array.shape[0], nvec, val=array)
  575 +
  576 +
309 577 if __name__ == '__main__':
310 578
311 579 # Demo the HarwellBoeingMatrix and RutherfordBoeingData classes
@@ -314,9 +582,7 @@ def readMatrix(self, fp, **kwargs):
314 582
315 583 import sys
316 584 fname = sys.argv[1]
317   - #M = HarwellBoeingMatrix(fname, patternOnly=False, readRhs=True)
318   - M = RutherfordBoeingData(fname, patternOnly=False)
319   - print 'Data of order (%-d,%-d) with %-d nonzeros' % (M.nrow,M.ncol,M.nnzero)
  585 + plot = False
320 586
321 587 numpy.set_printoptions(precision=8,
322 588 threshold=10,
@@ -324,26 +590,38 @@ def readMatrix(self, fp, **kwargs):
324 590 linewidth=80,
325 591 suppress=False)
326 592
  593 + #M = HarwellBoeingMatrix(fname, patternOnly=False, readRhs=True)
327 594 # Comment this out for Rutherford-Boeing data
328 595 #if M.readRhs:
329 596 # for i in range(M.nrhs):
330 597 # print M.rhs[:,i]
331 598
  599 + M = RutherfordBoeingData(fname, patternOnly=False)
  600 + print 'Data of order (%-d,%-d) with %-d nonzeros' % (M.nrow,M.ncol,M.nnzero)
  601 +
332 602 # Plot sparsity pattern
333   - try:
334   - import pylab
335   - except:
336   - sys.stderr.write('Pylab is required for the demo\n')
337   - sys.exit(1)
338   -
339   - (row,col) = M.find()
340   - fig = pylab.figure()
341   - ax = fig.gca()
342   - ax.plot(col, row, 'ks', markersize=1, linestyle='None')
343   - if M.issym: ax.plot(row, col, 'ks', markersize=1, linestyle='None')
344   - ax.set_xlim(xmin=-1, xmax=M.ncol)
345   - ax.set_ylim(ymin=M.nrow, ymax=-1)
346   - if M.nrow == M.ncol:
347   - ax.set_aspect('equal')
348   - pylab.title(M.title, fontsize='small')
349   - pylab.show()
  603 + if plot:
  604 + try:
  605 + import pylab
  606 + except:
  607 + sys.stderr.write('Pylab is required for the demo\n')
  608 + sys.exit(1)
  609 +
  610 + (row,col) = M.find()
  611 + fig = pylab.figure()
  612 + ax = fig.gca()
  613 + ax.plot(col, row, 'ks', markersize=1, linestyle='None')
  614 + if M.issym: ax.plot(row, col, 'ks', markersize=1, linestyle='None')
  615 + ax.set_xlim(xmin=-1, xmax=M.ncol)
  616 + ax.set_ylim(ymin=M.nrow, ymax=-1)
  617 + if M.nrow == M.ncol:
  618 + ax.set_aspect('equal')
  619 + pylab.title(M.title, fontsize='small')
  620 + pylab.show()
  621 +
  622 + # Write data back to file
  623 + (ip, ind, val) = M.data()
  624 + print val
  625 + #write_rb('newmat.rb', M.nrow, M.ncol, ip, ind, val, symmetric=M.issym)
  626 + #x = numpy.ones(10)
  627 + #write_aux_from_numpy('x.rb', x)

0 comments on commit c412015

Please sign in to comment.
Something went wrong with that request. Please try again.