# GNU General Public License # # Program Tops - a stack-based computing environment # Copyright (C) 1999-2005 Dale R. Williamson # # Author and copyright holder of op4_master.c, loadop4.c, saveop4.c, op4.c: # Al Danial # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software Foundation, # Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # This file is automatically generated from the tops file # admin/templates/op4_master.c. Changes applied to this # file will be lost. # #include "sparse.h" #include "op4.h" #include "mex.h" #define MAX(a,b) ((a) > (b)) ? (a) : (b) def get_file_suffix(file,ext): """ constants would like to have the constant arrays below defined in op4.h but they cause 'multiple definition' errors """ op4_store_str = [ ['d', 'n', 0], ['s', '1', 0], ['s', '2', 0] ]; op4_words_per_term = [-1, 1, 2, 2, 4]; Nm_per_term = [1, # numbers per term[real data type] = 1 2]; # numbers per term[complex data type] = 2 op4_type_str = [ ['R', 'S', 0], ['R', 'D', 0], ['C', 'S', 0], ['C', 'D', 0] ]; # internal functions def op4_filetype(filename): """ Returns 1 if a text file 2 if a binary file native endian 3 if a binary file flipped endian 0 if none of the above or file read error """ DEBUG = 0 int word, type, nRead #char *w #w = (char *) &word; fp = open(filename, "rb"); # open in binary mode nRead = fread(&word, BYTES_PER_WORD, 1, fp); if DEBUG: print "op4_filetype first word is %d chars:[%c][%c][%c][%c]\n" %(word, w[0], w[1], w[2], w[3]); if(nRead == 1): # in units of words if(word >= 538976288): # If this is a text file the smallest possible value for # the integer comprising the first four bytes is # 538976288 which corresponds to four blank spaces. if((ISBLANK( w[0]) and ISBLANK( w[1]) and # \s\s\s\s ISBLANK( w[2]) and ISBLANK( w[3])) or (ISBLANK( w[0]) and ISBLANK( w[1]) and # \s\s\s\d ISBLANK( w[2]) and isdigit((int) w[3])) or (ISBLANK( w[0]) and ISBLANK( w[1]) and # \s\s\d\d isdigit((int) w[2]) and isdigit((int) w[3])) or (ISBLANK( w[0]) and isdigit((int) w[1]) and # \s\d\d\d isdigit((int) w[2]) and isdigit((int) w[3])) or (isdigit((int) w[0]) and isdigit((int) w[1]) and # \d\d\d\d isdigit((int) w[2]) and isdigit((int) w[3]))): type = 1; else: type = 0; ### elif(word == 24): # Fortran binary record header, type = 2; # native endian, for 6 word block elif(word == 402653184): # Fortran binary record header, type = 3; # opposite endian, for 6 word block else: # the first word indicates the file is not an .op4 type = 0; ### else: # unable to read four bytes from the file type = 0; ### fp.close(); return type; def op4_scan( filename, # in n_mat, # out number of matrices name, # out matrix names storage, # out 0=dn; 1=sp1; 2=sp2 nRow, # out number of rows nCol, # out number of columns nStr, # out number of strings nNnz, # out number of nonzero terms type, # out 1=RS; 2=RD; 3=CS; 4=CD form, # out matrix form 6=symm, etc digits, # out size of mantissa offset # out byte offset to matrix ): """ Count the number and determine storage type of matrices in the text op4 file by scanning for matrix, column & string headers. Driver to text and binary versions. """ #int i, filetype, endian,result # = 0 causes segfault filetype = op4_filetype(filename); # 1 if a text file # 2 if a binary file native endian # 3 if a binary file flipped endian # 0 if none of the above or file error if filetype==1: result = op4_scan_t(filename , # in n_mat , # out number of matrices name , # out matrix names storage , # out 0=dn; 1=sp1; 2=sp2 nRow , # out number of rows nCol , # out number of columns nStr , # out number of strings nNnz , # out number of nonzero terms type , # out 1=RS; 2=RD; 3=CS; 4=CD form , # out matrix form 6=symm, etc digits , # out size of mantissa offset); # out byte offset to matrix elif filetype==2 or filetype==3: endian = filetype - 2; result = op4_scan_b(filename , # in endian , # in 0=native 1=flipped n_mat , # out number of matrices name , # out matrix names storage , # out 0=dn; 1=sp1; 2=sp2 nRow , # out number of rows nCol , # out number of columns nStr , # out number of strings nNnz , # out number of nonzero terms type , # out 1=RS; 2=RD; 3=CS; 4=CD form , # out matrix form 6=symm, etc offset); # out byte offset to matrix if endian: for i in range(n_mat): digits[i] = -1; else: for i in range(n_mat): digits[i] = 0; ### ### else: result = 0; ### return result; def op4_scan_t(filename , # in n_mat , # out number of matrices name, # out matrix names storage, # out 0=dn; 1=sp1; 2=sp2 nRow, # out number of rows nCol, # out number of columns nStr, # out number of strings nNnz, # out number of nonzero terms type, # out 1=RS; 2=RD; 3=CS; 4=CD form, # out matrix form 6=symm, etc digits, # out size of mantissa offset # out byte offset to matrix ): """ Count the number and determine storage type of matrices in the text op4 file by scanning for matrix, column & string headers. Text version. Keep track of byte location, matrix name and matrix dims. Here are some typical headers: 10 10 6 2EYE 1P,3E23.16 10 -10 6 2EYE 1P,3E23.16 123456789x123456789x123456789x123456789x123456789x123456789x |-> optional from here on The first is for dense and sparse form 1, and the second is for sparse form 2. The headers are distinguished from numeric data in that there is no decimal point in the third column--a quality all text based numeric data has, as in -3.937729118E+02-4.689545312E+07 3.766712501E+07-6.581960491E+06 0.000000000E+00 """ DEBUG = 0 #char line[OP4_TXT_LINE_SIZE + 1]; #int n_text_cols, text_col_width, mixed, nwords, start_row, nRead; #int location; #char junk1[12], junk2[12], junk3[12]; if DEBUG: print "op4_scan_t file=[%s]\n" %(filename) n_mat = 0; fp = open(filename, "r"); if DEBUG: print "op4_scan_t opened ok\n" while (not feof(fp)): location = fp.tell(); fgets(line, OP4_TXT_LINE_SIZE, fp); if(strlen(line) >= 8): # actually the op4 file is illegal if line length < 8 lineType = op4_line_type(line) if lineType==OP4_TEXT_NUMERIC: break elif lineType==OP4_TEXT_HEADER: # {{{2 if(n_mat and not storage[n_mat-1]): # previous matrix stored in dense scheme nNnz[n_mat-1] = nRow[n_mat-1]nCol[n_mat-1]; ### storage[n_mat] = 0; # dense unless proven otherwise nStr[n_mat] = 0; nNnz[n_mat] = 0; if(line[41] == 'P'): nRead = sscanf(line, "%8d%8d%8d%8d%8s%[ 1P,]%1d%1[E]%d%1[.]%d", &nCol[n_mat], &nRow[n_mat], &form[n_mat], &type[n_mat], name[n_mat], junk1, # 1P, &n_text_cols, junk2, # E &text_col_width, junk3, # . &digits[n_mat]); if(nRead != 11): fp.close() return 0; ### else: n_text_cols = 5 text_col_width = 15 digits[n_mat] = 9 nRead = sscanf(line, "%8d%8d%8d%8d%8s", &nCol[n_mat], &nRow[n_mat], &form[n_mat], &type[n_mat], name[n_mat]) if(nRead != 5): fp.close() return 0; ### ### nRow[n_mat] = fabs(nRow[n_mat]) offset[n_mat] = location ++(n_mat); elif lineType==OP4_TEXT_NEW_COLUMN: # {{{2 if DEBUG: print "new col: %s" %(line) elif lineType==OP4_TEXT_STRING_2: # {{{2 sscanf(line, "%8d%8d", &nwords, &start_row); --nwords; nNnz[n_mat-1] += nwords/op4_words_per_term[type[n_mat-1]] if DEBUG: print "s2: nwords=%d nNnz=%d\n" %(nwords, nNnz[n_mat-1]); storage[n_mat-1] = 2; ++nStr[ n_mat-1]; elif lineType==OP4_TEXT_STRING_1: if DEBUG: print "s1: %s" %(line) sscanf(line, "%8d", &mixed); nwords = (mixed/65536) - 1; start_row = mixed - 65536 * (nwords + 1); nNnz[ n_mat-1] += nwords/op4_words_per_term[type[n_mat-1]]; storage[n_mat-1] = 1; ++nStr[n_mat-1]; elif lineType==OP4_TEXT_ERROR: fp.close() return 0; ### ### ### ### offset[n_mat] = fp.tell(); # last byte in the file fp.close() if not storage[n_mat-1]: # last matrix stored in dense scheme nNnz[n_mat-1] = nRow[n_mat-1]nCol[n_mat-1]; ### return 1; def op4_scan_b(filename, # in endian , # in 0=native 1=flipped n_mat , # out number of matrices name, # out matrix names storage, # out 0=dn; 1=sp1; 2=sp2 nRow , # out number of rows nCol , # out number of columns nStr , # out number of strings nNnz , # out number of nonzero terms type , # out 1=RS; 2=RD; 3=CS; 4=CD form , # out matrix form 6=symm, etc offset # out byte offset to matrix ): """ Count the number and determine storage type of matrices in the text op4 file by scanning for matrix, column & string headers. Binary version. """ DEBUG = 0 #int record_length, is_header, result # = 0 causes segfault #char my_name[9], next_byte; #long location; if DEBUG: print "op4_scan_b file=[%s]\n" %(filename) n_mat = 0; fp = open(filename, "rb"); if DEBUG: print "op4_scan_b opened ok\n"; while (not feof(fp)): location = fp.tell(); if DEBUG: printf("op4_scan_b location A =%ld\n", location); # Determining EOF is tricky because positioning the file pointer # via Fortran record lengths ends up putting the FP at the very # last byte in the file... without triggering EOF. Have to # explicitly look ahead one byte to get feof() to register . next_byte = fgetc(fp); if(feof(fp)): # ok I'm really done with this file break; else # oops, am not at EOF, put the byte back ungetc(next_byte, fp); op4_is_mat_header_b(fp , endian , &record_length , &is_header , my_name , &storage[n_mat] , &nRow[n_mat] , &nCol[n_mat] , &nStr[n_mat] , &nNnz[n_mat] , &type[n_mat] , &form[n_mat] , &offset[n_mat]); if DEBUG: print "op4_scan_b got RL=%d\n" %(record_length) fseek(fp, record_length + 2*BYTES_PER_WORD, SEEK_CUR) # skip header if DEBUG: print "op4_scan_b loc=%ld is_header=%d\n" %(fp.tell(), is_header) if is_header: strncpy(name[n_mat], my_name, 9); if storage[n_mat]: # sparse 1 or 2; need to count strings result = op4_count_str_b(fp, endian, type[n_mat], storage[n_mat], nRow[n_mat], nCol[n_mat], &nStr[n_mat], &nNnz[n_mat]); if not result: fp.close() n_mat = 0 return result ### ### ++(n_mat); ### if DEBUG: location = fp.tell(); print "op4_scan_b location B =%ld nStr=%d\n" %(location, nStr[n_mat-1]); ### ### fseek(fp, 0, SEEK_END); # Set the file pointer to the end of file. offset[n_mat] = fp.tell(); # If this isn't done ftell() returns bogus. if DEBUG: print "op4_scan_b end byte position=%ld\n" %(offset[n_mat]) fp.close() if not storage[n_mat-1]: # last matrix stored in dense scheme nNnz[n_mat-1] = nRow[n_mat-1]nCol[n_mat-1]; ### return 1; def op4_is_mat_header_b(fp, # {{{1 endian , # in 0=native 1=flipped record_length, # out is_header, # out 1=is matrix hdr; 0=isn't my_name , # out matrix name storage , # out 0=dn; 1=sp1; 2=sp2 nRow , # out number of rows nCol , # out number of columns nStr , # out number of strings nNnz , # out number of nonzero terms type , # out 1=RS; 2=RD; 3=CS; 4=CD form , # out matrix form 6=symm, etc offset): # out byte offset to matrix """ If this Fortran record contains a matrix header, reads the appropriate entries. The file pointer is rewound to where it was when this routine was called. """ DEBUG = 0; #int is_ascii, nCol_, nRow_, Form_, Type_, rec_len, # column_id, start_row, n_words, header[7]; #long location; #char name_ptr; is_header = 0; location = fp.tell(); fread(record_length, BYTES_PER_WORD, 1, fp); if DEBUG: print "op4_is_mat_header_b RL=%d location=%ld\n" %(record_length, location) if endian: record_length = flip_bytes_int(record_length) if DEBUG: print "op4_is_mat_header_b RL=%d location=%ld\n" %(record_length, location) if(record_length == 24): # candidate for a matrix header fread(header, BYTES_PER_WORD, 7, fp); # + trailing rec marker # words header[4],header[5] make up the ASCII name--is it valid? name_ptr = header[4]; is_ascii = 1; if(op4_valid_name(name_ptr)): strncpy(my_name, name_ptr,9) else: strncpy(my_name, "unnamed",9) ### nCol_ = header[0]; nRow_ = header[1]; Form_ = header[2]; Type_ = header[3]; if(endian): nCol_ = flip_bytes_int(nCol_) nRow_ = flip_bytes_int(nRow_) Form_ = flip_bytes_int(Form_) Type_ = flip_bytes_int(Type_) ### if DEBUG: print "op4_scan_b name=[%s] " %(my_name) print " is_ascii=%d nCols=%d nRows=%d Form=%d Type=%d\n" %(is_ascii,nCol_,nRow_,Form_,Type_) ### if(is_ascii and nCol_ > 0 and nRow_ and Form_ < 50 and Type_ < 10): # we have a winner; this is a matrix header is_header = 1 nRow = fabs(nRow_) nCol = nCol_ offset = location type = Type_ form = Form_ nStr = 0 nNnz = 0 fread(&rec_len , BYTES_PER_WORD, 1, fp) fread(&column_id, BYTES_PER_WORD, 1, fp) fread(&start_row, BYTES_PER_WORD, 1, fp) fread(&n_words , BYTES_PER_WORD, 1, fp) if(endian): rec_len = flip_bytes_int(rec_len) column_id = flip_bytes_int(column_id) start_row = flip_bytes_int(start_row) n_words = flip_bytes_int(n_words) ### if DEBUG: print "op4_scan_b rl=%2d c=%2d r=%2d nw=%2d\n" %(rec_len, column_id, start_row, n_words) if(nRow_ < 0): storage = 2; # sparse 2 elif(not start_row or ((start_row == 1) and (n_words == 1))): # null storage = 1; # sparse 1 else: # dense storage = 0; ### ### ### fseek(fp, location, SEEK_SET); def op4_count_str_b(FILE fp , # in {{{1 int endian , # in int nType , # in int storage, # in int nRows , # in int nCols , # in int nStr , # out int nNnz): # out DEBUG = 0 #int col, row, n_words, length_in_words, record_length, encoded, word_count, end_of_data; #define ErrStrSize 200 #char error_msg[ErrStrSize - 1]; nStr = 0 col = 1 while (col <= nCols): if DEBUG: print "op4_count_str_b top of loop loc=%ld\n" %(fp.tell()) fread(&record_length, BYTES_PER_WORD, 1, fp); fread(&col, BYTES_PER_WORD, 1, fp); fread(&end_of_data, BYTES_PER_WORD, 1, fp); fread(&n_words, BYTES_PER_WORD, 1, fp); if endian: record_length = flip_bytes_int(record_length) col = flip_bytes_int(col) end_of_data = flip_bytes_int(end_of_data) n_words = flip_bytes_int(n_words) ### if DEBUG: printf("record_length =%d\n", record_length); printf("col =%d\n", col ); printf("end_of_data =%d\n", end_of_data ); printf("n_words =%d\n", n_words ); ### # make sure values aren't bogus {{{2 if((end_of_data != 0) and (end_of_data != 1)): nNnz = 0; snprintf(error_msg, ErrStrSize, "end_of_data=%d out of range, byte offset %ld\n", end_of_data, fp.tell()); raise Exception(error_msg); ### if((col < 1) or (col > (nCols + 2))): nNnz = 0; snprintf(error_msg, ErrStrSize, "col=%d out of range, byte offset %ld\n", col, fp.tell()); raise Exception(error_msg); ### if((record_length < 4*BYTES_PER_WORD) or (record_length > BYTES_PER_WORD * (5 + op4_words_per_term[nType]nRows))): # 5 = 3 column header ints + # 2 string header ints nNnz = 0; snprintf(error_msg, ErrStrSize, "RL=%d out of range (should be < %d), byte offset %ld\n", record_length, BYTES_PER_WORD*(5 + op4_words_per_term[nType]nRows), fp.tell()); raise Exception(error_msg); ### # 2}}} if DEBUG: printf("op4_count_str_b A RL=%d col=%d enddata=%d n_words=%d loc=%ld\n" %(record_length, col, end_of_data, n_words, fp.tell()); if end_of_data: fseek(fp, record_length - 2*BYTES_PER_WORD, SEEK_CUR); break; ### word_count = 0; while (word_count < n_words): if(storage == 1): fread(&encoded, BYTES_PER_WORD, 1, fp); word_count += 1; if endian: encoded = flip_bytes_int(encoded); ### length_in_words = (encoded/65536) -1; if DEBUG: print "op4_count_str_b B sp1 length_words=%d encoded=%d loc=%ld\n" %(length_in_words, encoded, fp.tell()) ### else: # storage == 2 fread(&length_in_words, BYTES_PER_WORD, 1, fp); fread(&row, BYTES_PER_WORD, 1, fp); # unused word_count += 2; if endian: length_in_words = flip_bytes_int(length_in_words); --length_in_words; if DEBUG: printf("op4_count_str_b C sp2 length_words=%d row=%d loc=%ld\n" %(length_in_words, row, fp.tell()); if(length_in_words < 0): snprintf(error_msg, ErrStrSize, "length_in_words=%d out of range, byte offset %ld\n", length_in_words, fp.tell()); raise Exception(error_msg); ### nNnz += length_in_words; ++(nStr); # skip numeric data fseek(fp, length_in_words*BYTES_PER_WORD, SEEK_CUR); word_count += length_in_words; if DEBUG: print "op4_count_str_b D word_count=%d loc=%ld\n" %(word_count, fp.tell()); fread(&record_length, BYTES_PER_WORD, 1, fp); if endian: record_length = flip_bytes_int(record_length); ### if DEBUG: print "op4_count_str_b E end RL=%d loc=%ld\n" %(record_length, fp.tell()) ### ### nNnz /= op4_words_per_term[nType]; return 1 # 1}}} def op4_line_type(line): # in {{{1 #int length, type; DEBUG = 0 length = len(line) if(length > 2 and line[2] == '.'): # most common case: numeric data type = OP4_TEXT_NUMERIC; elif(length >= 34): # header string must be at least 34 chars type = OP4_TEXT_HEADER; elif(length >= 24 and line[23] >= '0' and line[23] <= '9'): type = OP4_TEXT_NEW_COLUMN; elif(length >= 16 and line[15] >= '0' and line[15] <= '9'): type = OP4_TEXT_STRING_2; elif(length >= 8 and line[ 7] >= '0' and line[ 7] <= '9'): type = OP4_TEXT_STRING_1; else: if DEBUG: print "op4_line_type Error on [%s]\n" %(line) type = OP4_TEXT_ERROR if DEBUG: print "T=%d:[%s]\n", %(type,line) ### ### return type; # 1}}} def op4_read_col_t(fp , # {{{1 c_in , # in requested column to read nRow , # in # rows in matrix nCol , # in # columns in matrix fmt_str , # in eg "%23le%23le%23le" col_width , # in # characters in format str storage , # in 0=dense 1,2=sparse complx , # in 0=real 1=complex n_str , # in/out idx str_data[] (s_o)=1,2 S , # out string data (s_o)=1,2 N_index , # in/out idx N[] (s_o)=1,2 N # out numeric data DEBUG = 0 #double x[OP4_MAX_TEXT_COLS]; #int nRead, j, nNum_cols, still_reading, location, col, row, # nNum, got_this_col, n_nnz, mixed, nwords, start_row, # line_type, NPT, # n_read_this_str = 0; #char line[OP4_TXT_LINE_SIZE + 1]; still_reading = 1; got_this_col = 0; n_nnz = 0; # count of numeric values loaded NPT = Nm_per_term[complx]; # numbers per term; real=1 complex=2 if not storage: # zero out dense column for j in range(NPTnRow): N[j] = 0.0; ### while (still_reading): location = fp.tell(); fgets(line, OP4_TXT_LINE_SIZE, fp); line_type = op4_line_type(line); if DEBUG: print "op4_read_col_t line_type=%d\n" %(line_type) if(strlen(line) >= 8): switch (line_type) case OP4_TEXT_NUMERIC : # {{{2 nNum_cols = (int) (len(line) / col_width) nRead = sscanf(line, fmt_str, &x[0], &x[1], &x[2], &x[3], &x[4], &x[5], &x[6], &x[7], &x[8]); if storage: # sparse if DEBUG: print "op4_read_col_t inserting %d terms at row=%d\n" %(nRead, S[n_str-1].start_row) n_read_this_str += nRead; S[n_str-1].len = n_read_this_str/NPT; for j in range(nRead): N[(N_index)++] = x[j]; else: # dense if DEBUG: print "op4_read_col_t inserting %d terms at row=%d\n" %(nRead, row); for j in range(nRead): if DEBUG: print " N[row=%d + N_index=%d]=%e\n" %(row, (N_index), x[j]) N[row + j] = x[j]; ### row += nRead; ### n_nnz += nRead; if DEBUG: print "op4_read_col_t nR=%d n_nnz now=%d\n" %(nRead, n_nnz); case OP4_TEXT_HEADER : # {{{2 raise Exception("op4_read_col_t new matrix header while reading column"); case OP4_TEXT_NEW_COLUMN: # {{{2 if DEBUG: print "op4_read_col_t start new column (got_this_col=%d) loc=%ld\n" %(got_this_col, fp.tell()); if(got_this_col): # have read next column's header; back up & exit fseek(fp, location, SEEK_SET); still_reading = 0; else: nRead = sscanf(line, "%8d%8d%8d", &col, &row, &nNum); if DEBUG: print " requested col=%d this col=%d this row=%d\n", c_in, col, row) if(c_in != col): # this is not the column of interest fseek(fp, location, SEEK_SET); return n_nnz; ### --col; --row; # internal indexing is 0 based if not storage: # dense row *= NPT if DEBUG: print "\n new column %2d row %2d complx=%d\n" %(col, row, complx); got_this_col = 1; if((col > nCol) or (nRead != 3)): # have reached the spurious end of matrix entry # or have encountered bad data return n_nnz; ### ### case OP4_TEXT_STRING_1 : case OP4_TEXT_STRING_2 : # {{{2 n_read_this_str = 0; if(line_type == OP4_TEXT_STRING_1): sscanf(line, "%8d", &mixed); nwords = (mixed/65536) - 1; start_row = mixed - 65536 * (nwords + 1); else: # type 2 sscanf(line, "%8d%8d", &nwords, &start_row); --nwords; ### S[n_str].start_row = start_row - 1; S[n_str].len = 0; S[n_str].N_idx = N_index; ++(n_str); case OP4_TEXT_ERROR: # {{{2 raise Exception("op4_read_col_t error reading a column "); ### ### ### return n_nnz; # 1}}} def op4_read_col_b(fp , # {{{1 endian , # in 0=native 1=flipped c_in , # in requested column to read nRow , # in # rows in matrix nCol , # in # columns in matrix nType , # in 1=RS 2=RD 3=CS 4=CD storage , # in 0=dn; 1=sp1; 2=sp2 n_str , # in/out idx str_data[] (s_o)=1,2 S , # out string data (s_o)=1,2 N_index , # in/out idx N[] (s_o)=1,2 N # out numeric data ): DEBUG = 0 #int complx, NPT, WPN, BPN, i, record_length, col, start_row, # n_numbers, got_this_col, n_nnz, is_header, unused, mixed, # n_w_this_col, n_w_this_str; #long l_unused; #char unused_name[9]; #float xs; if DEBUG: print "\nop4_read_col_b 1 location=%ld endian=%d want col %d\n" %(fp.tell(), endian, c_in); complx = nType > 2 ? 1 : 0; NPT = Nm_per_term[complx]; # numbers per term; real=1 complex=2 WPN = nType % 2 ? 1 : 2; # words per number; float=1 double=2 BPN = BYTES_PER_WORD * WPN; # bytes per number got_this_col = 0; n_nnz = 0; # count of numeric values loaded op4_is_mat_header_b(fp , # in endian , # in &record_length , # out &is_header , # out unused_name , # out &unused , # out &unused , # out &unused , # out &unused , # out &unused , # out &unused , # out &unused , # out &l_unused); # out if DEBUG: print "op4_read_col_b 2 RL=%d location=%ld\n" %(record_length, fp.tell()) if(is_header): # yes? skip it fseek(fp, record_length+2*BYTES_PER_WORD, SEEK_CUR); if DEBUG: print "op4_read_col_b 3 RL=%d location=%ld\n" %(record_length, fp.tell()) if not storage: # zero out dense column for i in range(NPTnRow): N[i] = 0.0; ### ### fread(&record_length, BYTES_PER_WORD, 1, fp); fread(&col, BYTES_PER_WORD, 1, fp); fread(&start_row, BYTES_PER_WORD, 1, fp); fread(&n_w_this_col, BYTES_PER_WORD, 1, fp); if endian: record_length = flip_bytes_int(record_length); col = flip_bytes_int(col); start_row = flip_bytes_int(start_row); n_w_this_col = flip_bytes_int(n_w_this_col); ### if DEBUG: print "op4_read_col_b: requested col=%d RL=%d got col=%d row=%d nW=%d\n" %(c_in, record_length, col, start_row, n_w_this_col) n_numbers = n_w_this_col/WPN; if(c_in > col): # requested column ahead of me; skip record fseek(fp, record_length-2*BYTES_PER_WORD, SEEK_CUR); if DEBUG: print "op4_read_col_b: column mismatch, advance to location %ld & exit\n" %(fp.tell()); return 0; elif(c_in < col): # requested column behind me; back up fseek(fp,-4*BYTES_PER_WORD, SEEK_CUR); if DEBUG: print "op4_read_col_b: column mismatch, back up to location %ld & exit\n" %(fp.tell()) return 0; ### if DEBUG: print "op4_read_col_b: going to read col=%d\n" %(col) if not storage: # dense; read entire column if(WPN == 2): # double precision: can do bulk read if DEBUG: print "op4_read_col_b: double prec\n" fread(&N[NPT*(start_row-1)], BPN, n_numbers, fp); if endian: # opposite endian, need to flip terms for i in range(n_numbers): N[NPT*(start_row-1) + i] = flip_bytes_double(N[NPT*(start_row-1) + i]); ### ### ### elif not endian: # single precision native endian if DEBUG: print "op4_read_col_b: single prec native endian\n" for i in range(n_numbers): fread(&xs, BPN, 1, fp); N[NPT*(start_row-1) + i] = xs; ### else: # single precision opposite endian if DEBUG: print "op4_read_col_b: single prec opposite endian\n" for i in range(n_numbers): fread(&xs, BPN, 1, fp); N[NPT*(start_row-1) + i] = flip_bytes_float(xs); ### ### if DEBUG: for i in range(n_numbers): print "op4_read_col_b: N[%d]=%le\n" %(NPT*(start_row-1) + i,N[NPT*(start_row-1) + i]) ### ### fseek(fp, BYTES_PER_WORD, SEEK_CUR); # skip trailing RL if DEBUG: print "op4_read_col_b: finished numeric data location=%ld\n" %(fp.tell()); else: # sparse; loop over all strings record_length -= 2*BYTES_PER_WORD; # don't count start, end rec mark while (record_length >= 8): if(storage == 1): # sparse 1 fread(&mixed, BYTES_PER_WORD, 1, fp); if endian: mixed = flip_bytes_int(mixed); n_w_this_str = (mixed/65536) - 1; start_row = mixed - 65536 * (n_w_this_str + 1); record_length -= 1*BYTES_PER_WORD; else: # sparse 2 fread(&n_w_this_str, BYTES_PER_WORD, 1, fp); fread(&start_row, BYTES_PER_WORD, 1, fp); record_length -= 2*BYTES_PER_WORD; if endian: n_w_this_str = flip_bytes_int(n_w_this_str); start_row = flip_bytes_int(start_row); ### --n_w_this_str; ### n_numbers = n_w_this_str/WPN; if DEBUG: print "\nop4_read_col_b top of str: start_row=%d n_w_this_str=%d n_numbers=%d N_index=%d\n" %(start_row, n_w_this_str, n_numbers, N_index) S[n_str].start_row = start_row - 1; S[n_str].len = n_numbers/(complx + 1); S[n_str].N_idx = N_index; ++(n_str); if(WPN == 2): # double precision: can do bulk read if DEBUG: print "op4_read_col_b sp: double prec\n" fread(&N[N_index], BPN, n_numbers, fp); if endian: # opposite endian, need to flip terms for i in range(n_numbers): N[N_index + i] = flip_bytes_double(N[N_index + i]); ### ### elif not endian: # single precision native endian if DEBUG: print "op4_read_col_b sp: single prec native endian\n") for i in range(n_numbers): fread(&xs, BPN, 1, fp); N[N_index + i] = xs; ### else: # single precision opposite endian if DEBUG: print "op4_read_col_b sp: single prec opposite endian\n" for i in range(n_numbers): fread(&xs, BPN, 1, fp); N[N_index + i] = flip_bytes_float(xs); ### ### record_length -= BYTES_PER_WORD*n_w_this_str; if DEBUG: print "op4_read_col_b sp end string RL=%d loc=%ld\n" %(record_length, fp.tell()); for i in range(n_numbers): print "op4_read_col_b N[%3d] = %e\n" %(N_index+i, N[N_index + i]) ### ### N_index += n_numbers n_nnz += n_numbers ### fread(&record_length, BYTES_PER_WORD, 1, fp) ### return n_nnz; FILE* op4_open_r(const char filename, long offset): # in FILE fp; # Open an op4 file for reading, optionally positioning the file # handle to a location past the start of the file. # fp = fopen(filename, "r"); if(fp == NULL) return NULL; if(offset) fseek(fp, offset, SEEK_SET); return fp; def int flip_bytes_int(int x): #int y; #char *c_x, *c_y; c_x = (char *) &x; c_y = (char *) &y; c_y[0] = c_x[3]; c_y[1] = c_x[2]; c_y[2] = c_x[1]; c_y[3] = c_x[0]; return y; def float flip_bytes_float(float x): """ Looks identical to flip_bytes_int() and it practice it might be possible to use this function and flip_bytes_int() interchangably. The combination of gcc and x86 hardware may cause problems though since floats might be assigned extra bytes internally. """ float y; char *c_x, *c_y; c_x = (char *) &x; c_y = (char *) &y; c_y[0] = c_x[3]; c_y[1] = c_x[2]; c_y[2] = c_x[1]; c_y[3] = c_x[0]; return y; def double flip_bytes_double(double x): double y; char *c_x, *c_y; c_x = (char *) &x; c_y = (char *) &y; c_y[0] = c_x[7]; c_y[1] = c_x[6]; c_y[2] = c_x[5]; c_y[3] = c_x[4]; c_y[4] = c_x[3]; c_y[5] = c_x[2]; c_y[6] = c_x[1]; c_y[7] = c_x[0]; return y; def op4_valid_name(name): """ Returns: 1 if the name passed in is a valid Nastran op4 matrix name -1 if the name is null 0 if the name is not valid @todo this check is too conservative; Nastran allows names like 'a+' """ DEBUG = 0 #int i, len, fixed_len; if DEBUG: print "op4_valid_name with [%s]\n" %(name) length = len(name) # strip trailing blanks fixed_len = len; name = name.strip() len = len(name) if not length: # null string? return -1; elif(length > 8): # too long? return 0; elif(not isalpha((int) name[0])): # first char must be [a-zA-Z] return 0; ### for i in range(len): # this test is too conservative; Nastran allows names like 'a+' if(not (ISUNERDSCORE( name[i]) or isalpha((int) name[i]) or isdigit((int) name[i]))): return 0; ### ### return 1; # words def mexFunction(nlhs, plhs[], int nrhs, prhs[]): DEBUG = 0 filename[80], line[OP4_TXT_LINE_SIZE+1], name[Max_Matrices_Per_OP4][9], fmt_str[OP4_FMT_LINE_SIZE], fmt_one[OP4_FMT_STR_SIZE]; int c, i, j, s_ptr, n_ptr, n_nnz, n_mat_in_file, complx, rows, NPT, Nmat, Nskip, nNum_cols, filetype, col_width = 0, load_all = 0, unused = 0, storage[Max_Matrices_Per_OP4], nRow[Max_Matrices_Per_OP4], nCol[Max_Matrices_Per_OP4], nStr[Max_Matrices_Per_OP4], nNnz[Max_Matrices_Per_OP4], type[Max_Matrices_Per_OP4], form[Max_Matrices_Per_OP4], digits[Max_Matrices_Per_OP4]; long offset[Max_Matrices_Per_OP4]; #double nskip_ptr; double Mr[Max_Matrices_Per_OP4], Mi[Max_Matrices_Per_OP4]; int ir[Max_Matrices_Per_OP4], jc[Max_Matrices_Per_OP4]; #int r, s, ir_index, n_str_one_col, # ptr_H, ptr_S_start, ptr_N_start, ptr_S, ptr_N, size print_header = 1; # 1=print matrix info after loading; 0=don't SparseMatrix m; D = 0.0; if((nrhs < 0) or (nrhs > 2)): raise Exception("Can only have one or two RHS arguments.\n") if(not mxIsChar(prhs[0])): raise Exception("First RHS argument must be a text filename.\n") mxGetString(prhs[0], filename, 80) Nmat = nlhs if((Nmat < 1) or (Nmat >= Max_Matrices_Per_OP4) ): raise Exception("Number of LHS arguments must be >= 1 and <= 50.\n") ### if(nrhs > 1): # get # of matrices to skip over nskip_ptr = mxGetPr(prhs[1]) Nskip = (int) nskip_ptr if(Nskip < 0): raise Exception("Second RHS argument must be an integer >= 0.\n") ### else: Nskip = 0 ### if DEBUG: print "loadop4 nskip=%d nmat=%d name=[%s]\n" %(Nskip, Nmat, filename); # Scan the file for number of matrices and headers for each matrix. if(not op4_scan(filename, # in &n_mat_in_file, # out number of matrices name, # out matrix names storage, # out 0=dn; 1=sp1; 2=sp2 nRow, # out number of rows nCol, # out number of columns nStr, # out number of strings nNnz, # out number of nonzero terms type, # out 1=RS; 2=RD; 3=CS; 4=CD form, # out matrix form 6=symm, etc digits, # out size of mantissa offset)): # out byte offset to matrix snprintf(line, OP4_TXT_LINE_SIZE, " file error for '%s'", filename); raise Exception(line); ### if DEBUG: print "loadop4 n_mat_in_file = %d\n" %(n_mat_in_file) if(n_mat_in_file < 1 ): snprintf(line, OP4_TXT_LINE_SIZE, " no matrices in '%s'", filename); raise Exception(line); return if(Nmat == -1): # wants every matrix in the file (beyond nSkip) load_all = 1 Nmat = n_mat_in_file - Nskip if((Nmat < 0) or (Nskip + Nmat) > n_mat_in_file): # Requested more matrices than exist in file. snprintf(line, OP4_TXT_LINE_SIZE, " only %d matrices in '%s'",n_mat_in_file, filename); raise Exception(line); ### filetype = op4_filetype(filename); fp = op4_open_r(filename, 0); if(fp == NULL): raise Exception("loadop4: unable to open op4 file at offset\n"); ### for i in range(Nskip,Nskip+Nmat): if DEBUG: print "%2d %2d x %2d %2s %2s nStr=%4d nNnz=%4d f=%2d %-8s D=%2d O=%12ld L=%ld\n" %(i+1, nRow[i], nCol[i], op4_type_str[type[i]-1], op4_store_str[storage[i]], nStr[i], nNnz[i], form[i], name[i], digits[i], offset[i], offset[i+1] - offset[i]) # allocate memory, place stack item to receive op4 matrix {{{2 complx = type[i] > 2 ? 1 : 0; NPT = Nm_per_term[complx]; if storage[i]: # sparse # matlab sparse {{{3 # Declare a tops native sparse matrix having only 1 column. # Assume worst case values for # strings, # non-zeros. Have # to declare a double to force alignment to multiple of 8. # This code duplicates code in the following functions: # - src/sparce.c:sp_data_offsets() # - src/sparce.c:sp_set_header() # - src/sparce.c:sparse_overlay() n_str_one_col = nRow[i]/2 + 1; ptr_H = SPARSE_MAGIC_SIZE; ptr_S_start = ptr_H + SPARSE_HDR_SIZE; ptr_N_start = ptr_S_start + (1 + 1); # nCOl = 1 ptr_S = ptr_N_start + (1 + 1); # nCOl = 1 ptr_N = ptr_S + n_str_one_col*sizeof(str_t)/ sizeof(int); if(ptr_N % 2): ++(ptr_N) size = (ptr_N)*sizeof(int) + NPTnRow[i]*sizeof(double) if((D = (double *) malloc(size)) == NULL): raise Exception("loadop4: insufficient memory for matrix column"); ### m.data = (char *) D; m.magic = (int *) m.data; m.H = &m.magic[SPARSE_MAGIC_SIZE]; m.S_start = &m.magic[ptr_S_start]; m.N_start = &m.magic[ptr_N_start]; m.S = (str_t *) &m.magic[ptr_S ]; m.N = (double *) &m.magic[ptr_N ]; m.H[ROWS] = nRow[i]; m.H[COLS] = 1; m.H[COMPLX] = complx; m.H[n_STR] = nRow[i]/2 + 1; m.H[n_NONZ] = nRow[i]*NPT; m.H[DATA_SIZE]= size; # then create the matlab variables if complx: # sparse complex plhs[i-Nskip] = mxCreateSparse(nRow[i], nCol[i], nNnz[i], mxCOMPLEX); Mr[i-Nskip] = mxGetPr(plhs[i-Nskip]); Mi[i-Nskip] = mxGetPi(plhs[i-Nskip]); else: # sparse real plhs[i-Nskip] = mxCreateSparse(nRow[i], nCol[i], nNnz[i], mxREAL); Mr[i-Nskip] = mxGetPr(plhs[i-Nskip]); ### ir[i-Nskip] = mxGetIr(plhs[i-Nskip]); jc[i-Nskip] = mxGetJc(plhs[i-Nskip]); # 3}}} else: # dense if complx: # dense complex plhs[i-Nskip] = mxCreateDoubleMatrix(nRow[i], nCol[i], mxCOMPLEX); Mr[i-Nskip] = mxGetPr(plhs[i-Nskip]); Mi[i-Nskip] = mxGetPi(plhs[i-Nskip]); for j in range(nRow[i]*nCol[i]): Mi[i-Nskip][j] = 0.0; else: # dense real plhs[i-Nskip] = mxCreateDoubleMatrix(nRow[i], nCol[i], mxREAL); Mr[i-Nskip] = mxGetPr(plhs[i-Nskip]); ### for j in range(nRow[i]*nCol[i]): Mr[i-Nskip][j] = 0.0; ### ### # 2}}} fseek(fp, offset[i], SEEK_SET); # jump to the ith matrix if(filetype == 1): # text fgets(line, OP4_TXT_LINE_SIZE, fp); # skip the header line # create the scanf format string, eg col_width = digits[i] + 7; strncpy( fmt_str, "", OP4_FMT_STR_SIZE); snprintf(fmt_one, OP4_FMT_STR_SIZE, "%%%dle", col_width); nNum_cols = (int) (80 / col_width); for c in range(nNum_cols): strncat(fmt_str, fmt_one, OP4_FMT_STR_SIZE); if DEBUG: print "one=[%s] str=[%s]\n" %(fmt_one, fmt_str) if storage[i]: # sparse jc[i-Nskip][0] = 0; ir_index = 0; else: # dense if((D = (double *) malloc(nRow[i]*NPT*sizeof(double))) == NULL): raise Exception("loadop4: insufficient memory for matrix column"); ### ### for c in range(1,nCols[i]): # op4 columns start at 1 if DEBUG: print "loadop4->column %2d n_ptr=%d\n" %(c, n_ptr) # read column c {{{2 if storage[i]: # sparse s_ptr = 0; n_ptr = 0; m.S_start[0] = 0; m.N_start[0] = 0; m.S[0].start_row = 0; m.S[0].len = 0; m.S[0].N_idx = 0; ### if(filetype == 1): # text if storage[i]: # load as sparse matrix n_nnz = op4_read_col_t(fp, c, nRow[i], nCol[i], fmt_str, col_width, storage[i], # in 0=dn 1,2=sp1,2 3=ccr complx , # in 0=real 1=complex &s_ptr , # in/out index m.S (if sp1,2) m.S , # out string data (if sp1,2) &n_ptr , # in/out index m.N m.N); # out numeric data else: # load as dense matrix n_nnz = op4_read_col_t(fp, c, nRow[i], nCol[i], fmt_str, col_width, storage[i], # in 0=dn 1,2=sp1,2 3=ccr complx , # in 0=real 1=complex unused , # in/out index m.S (if sp1,2) (str_t *) unused , # out string data (if sp1,2) (int *) unused , # in/out index m.N (sp 1,2) D); # out numeric data ### else: # binary if storage[i]: # load as sparse matrix n_nnz = op4_read_col_b(fp, filetype - 2, c, nRow[i], nCol[i], type[i], storage[i], # in 0=dn 1,2=sp1,2 &s_ptr , # in/out index m.S (if sp1,2) m.S , # out string data (if sp1,2) &n_ptr , # in/out index m.N m.N); # out numeric data else: # load as dense matrix n_nnz = op4_read_col_b(fp, filetype - 2, c, nRow[i], nCol[i], type[i], storage[i], # in 0=dn 1,2=sp1,2 unused , # in/out index m.S (if sp1,2) (str_t *) unused , # out string data (if sp1,2) (int *) unused , # in/out index m.N (sp 1,2) D ); # out numeric data ### ### # 2}}} # copy column c from internal to matlab variable if storage[i]: # sparse for s in range(s_ptr): for j in range(m.S[s].len): ir[i-Nskip][ir_index] = m.S[s].start_row + j; Mr[i-Nskip][ir_index] = m.N[m.S[s].N_idx + j*NPT]; if complx: Mi[i-Nskip][ir_index] = m.N[m.S[s].N_idx + j*NPT+1]; ++ir_index; ### ### jc[i-Nskip][c] = ir_index; else: # dense for r in range(nRow[i]): Mr[i-Nskip][(c-1)nRow[i] + r] = D[r*NPT]; if complx: Mi[i-Nskip][(c-1)nRow[i] + r] = D[r*NPT + 1]; ### ### if DEBUG: print "loadop4<-column %2d read %d values\n\n" %(c, n_nnz) ### ### end loop over columns of matrix i if(not storage[i]) # dense matrix free(D); if(print_header): print "loadop4: %2d x %2d %2s %2s nStr=%4d nNnz=%4d f=%2d %-8s dg=%1d of %1d\n" %(nRow[i], nCol[i], op4_type_str[type[i]-1],op4_store_str[storage[i]], nStr[i], nNnz[i],form[i],name[i], digits[i], i+1, Nmat) ### ### # end loop over matrices fp.close() if DEBUG: print "end of loadop4\n" return; ###