diff --git a/src/config/makefile.h b/src/config/makefile.h index 1669144..522dd7b 100644 --- a/src/config/makefile.h +++ b/src/config/makefile.h @@ -161,15 +161,7 @@ ifeq ($(EXPLICITF),TRUE) # with CPP to get .f files # .SUFFIXES: -.SUFFIXES: .o .s .F .f .c - -.F.o: - $(MAKE) $*.f - $(FC) -c $(FFLAGS) $*.f - /bin/rm -f $*.f - -.F.f: - $(CPP) $(INCLUDES) $(DEFINES) < $*.F | sed '/^#/D' > $*.f +.SUFFIXES: .o .s .c .c.o: $(CC) $(CFLAGS) -c $*.c diff --git a/src/inp/GNUmakefile b/src/inp/GNUmakefile index 94c090a..17383c8 100644 --- a/src/inp/GNUmakefile +++ b/src/inp/GNUmakefile @@ -3,7 +3,7 @@ OBJ = inp.o LIBRARY = libinp.a LIB_TARGETS = test.o test - HEADERS = inp.fh + HEADERS = inp.h include ../config/makefile.h diff --git a/src/inp/inp.c b/src/inp/inp.c new file mode 100644 index 0000000..9cb1eb7 --- /dev/null +++ b/src/inp/inp.c @@ -0,0 +1,630 @@ +#include +#include + +#include "inpP.h" + +int iread = 5; +int iwrite = 6; +int jrec = -1; +int jump = 0; +bool oswit = false; +int nerr = 999; +int nline = 0; +int noline = 0; +int ierrpos = -1; +char errmsg[2] = " "; +int input_line = 0; +char xblnk[2] = " "; +char xtab[2] = "\t"; +char xsplit[2] = ";"; +char xcomm[2] = "#"; +char xback[2] = "\\"; +char xquote[2] = "\""; + + +void inp_init(int ir, int iw) { + iread = ir; + iwrite = iw; + jrec = -1; + jump = 0; + oswit = false; + nerr = 999; + nline = 0; + noline = 0; + ierrpos = -1; + errmsg[0] = ' '; + input_line = 0; +} + +int inp_n_field() { + + // return no. of fields in the input line ... 0 = EOF + return jump; +} + +int inp_cur_field() { + + // return no. of fields processed so far (0,...,inp_n_field()) + return jrec; +} + +void inp_set_field(int ivalue) { + + // set field to be read next (ivalue=0,...,inp_n_field()) + if (ivalue < 0 || ivalue > inp_n_field()) { + errquit('inp_set_field: stupid field value',ivalue); + } + jrec = ivalue; + + ierrpos = -1; + strcpy(errmsg, " "); + +} + +bool inp_line(char *z) { + + /* + set the variable z to be as much of the current input line + that it can hold + */ + + bool ret_val; + if (jump > 0) { + strcpy(z, ia); + ret_val = true; + ierrpos = -1; + strcpy(errmsg, " "); + } else { + strcpy(errmsg, "no input line available"); + ierrpos = -1; + ret_val = false; + } + + return ret_val; +} + +bool ois_ws(char xtest) { + return strcmp(xtest, xblnk) == 0 || strcmp(xtest, xtab) == 0; +} + +bool inp_read() { +/* + this routine reads a data card and scans it for non - space fields + the number of fields is stored in jump, the starting point of a + field in istrt(i) and the number of characters in that field + in inumb(i). +*/ + int ncol[MAX_FIELD], lenja, i, iwidth, j, jwidth, k, mark, nbegin, nfini; + bool ios_ws, ret_val; + char tmp[MAX_WIDTH], xprev; + + nline++; + if (nline <= noline) { + goto L150; + } + + if (oswit) { + ierrpos = -1; + strcpy(errmsg, "unexpected end of data file"); + jump = 0; + jrec = 0; + return false; + } else { + ierrpos = -1; + strcpy(errmsg, " "); + } + + // read next physical input line + lenja = 0; +L100: + scanf("%s", ja + lenja); + input_line++; + lenja = strlen(ja); + + // Check for . * eof at beginning of line to indicate EOF + if (lenja == 1 && (ja[0] == '.' || ja[0] == '*')) { + goto L300; + } + if (lenja == 3 && strcmp(ja, "eof") == 0) { + goto L300; + } + + // handle blank lines and concatenation using backslash + if (lenja == 0) { + goto L100; + } else { + if (ja[lenja - 1] == xback) { + ja[lenja - 1] = xblnk; + goto L100; + } + } + jwidth = strlen(ja); + + // handle comments from # to eol ... allow for backslash quoting + xprev = xblnk; + for (i = 0; i < jwidth; i++) { + if (ja[i] == xcomm && xprev != xback) { + lenja = strlen(ja); + printf("\n comment :- %s", ja + i + 1); + memset(ja + i, xblnk, MAX_WIDTH - i); + break; + } else if (ja[i] == xcomm && xprev == xback) { + strcpy(tmp, ja); + memmove(ja + i - 1, tmp + i, MAX_WIDTH - i); + xprev = xblnk; + continue; + } + strcpy(xprev, ja[i]); + } + + /* + figure out where ; splits physical line into multiple logical lines + again handling quoted backslash + */ + k = jwidth; + mark = 0; + xprev = xblnk; + for (i = 0; i < jwidth; i++) { + if (ja[i] == xsplit && xprev != xback) { + mark = mark + 1; + ncol[mark] = i; + } else if (ja[i] == xsplit && xprev == xback) { + strcpy(tmp, ja); + memmove(ja + i - 1, tmp + i, MAX_WIDTH - i); + xprev = xblnk; + continue; + } + xprev = ja[i]; + } + + noline = 1; + if (mark == 0) { + nstart[noline] = 1; + nend[noline] = jwidth; + } else { + i = ncol[mark] + 1; + if (i <= jwidth) { + for (j = i; j < jwidth; j++) { + if (!ois_ws(ja[j])) { + goto L170; + } + } + } + k = ncol[mark] - 1; + mark = mark - 1; +L170: + noline = mark + 1; + nstart[1] = 1; + for (i = 1; i <= mark; i++) { + j = ncol[i]; + nend[i] = j - 1; + nstart[i + 1] = j + 1; + } + nend[noline] = k; + } + nline = 1; + + // Start processing next logical input line (put into ia(1:iwidth)) +L150: + jump = 0; + jrec = 0; + nbegin = nstart[nline]; + nfini = nend[nline]; + iwidth = nfini - nbegin + 1; + memset(ia, xblnk, MAX_WIDTH); + memcpy(ia, ja + nbegin - 1, iwidth); + /* + partition input line into strings inside double quotes or + white space separated fields + */ + i = 1; +L151: + for (j = i; j <= iwidth; j++) { + if (!ois_ws(ia[j])) { + goto L152; + } + } + goto L155; +L152: + i = j; + jump++; + istrt[jump] = i; + if (ia[i] == xquote) { + for (j = i + 1; j <= iwidth; j++) { + if (ia[j] == xquote && ia[j - 1] != xback) { + break; + } else if (ia[j] == xquote && ia[j - 1] == xback) { + strcpy(tmp, ia); + memmove(ia + j - 1, tmp + j, MAX_WIDTH - j); + continue; + } + } + ierrpos = j; + strcpy(errmsg, "no terminating quote for string"); + return false; + } else { + for (j = i + 1; j <= iwidth; j++) { + if (ois_ws(ia[j])) { + break; + } + } + j--; + } + + inumb[jump] = j - istrt[jump] + 1; + i = j + 1; + goto L151; +L155: + if (jump > 0) { + iwidth = istrt[jump] + inumb[jump] - 1; + } + return true; +L300: + oswit = true; + ierrpos = -1; + strcpy(errmsg, "unexpected end of data file"); + jump = 0; + jrec = 0; + return false; +} + +bool inp_eof() { + return oswit; +} + +void inp_errout() { + char xpt, xstp; + xpt = '*'; + xstp = '.'; + int length, i; +/* + If an error has occured print out the error message + and the position in the current input line +*/ + if (strcmp(errmsg, " ") != 0) { + length = strlen(errmsg); + printf("input error at line %d: %s\n", input_line, errmsg); + jrec = -1; + printf("%s\n", ia); + if (ierrpos > 0) { + for (i = 0; i < ierrpos; i++) { + tmp[i] = xstp; + } + tmp[ierrpos] = xpt; + printf("%s\n", tmp); + } + } +} + +void inp_outrec() { + // Write out the current input line + printf("%s\n", ia); +} + +#include +#include + +bool inp_a(char* a) { +/* + Return field as character string, minus any enclosing quotes + with an error if it does not fit +*/ + int i1, i2, length; + bool ret_val; + + ierrpos = -1; + strcpy(errmsg, " "); + + if (jrec >= jump) { + a = xblnk; + ierrpos = 0; + strcpy(errmsg, "at end of line looking for character string"); + return false; + } + + i1 = istrt[jrec+1]; + i2 = istrt[jrec+1] + inumb[jrec+1] - 1; + + if (ia[i1] == xquote && ia[i2] == xquote) { + i1 = i1 + 1; + length = inumb[jrec+1] - 2; + } else { + length = inumb[jrec+1]; + } + + if (strlen(a) < length) { + a = xblnk; + ierrpos = 0; + strcpy(errmsg, "inp_a: string is too large for argument"); + return false; + } else { + jrec = jrec + 1; + strncpy(a, ia + i1, length); + return true; + } +} +logical function inp_a_trunc(a) + implicit none + + character*1 xblnk, xquote + character*(*) a +c +c Return field as character string, minus any enclosing quotes +c quietly truncating if it does not fit +c + ierrpos = -1 + errmsg = ' ' + if(jrec .ge. jump) then + a = xblnk + inp_a_trunc = .false. + ierrpos = 0 + errmsg = 'at end of line looking for character string' + return + endif + i1 = istrt(jrec+1) + i2 = istrt(jrec+1)+inumb(jrec+1)-1 + if (ia(i1:i1).eq.xquote .and. ia(i2:i2).eq.xquote) then + i1 = i1+1 + length = inumb(jrec+1)-2 + else + length = inumb(jrec+1) + endif + jrec = jrec + 1 + a = ia(i1:i1+length-1) + inp_a_trunc = .true. + return +c + end +logical function inp_f (buf) + implicit none + + double precision ten, buf + character*1 xchar(15) + data xchar /'0','1','2','3','4','5','6','7','8','9' + 1 ,'+','-','.','e','d'/ + data ten/10.0d0/ +c + ierrpos = -1 + errmsg = ' ' + buf=0.0d0 + if (jrec.ge.jump) then + inp_f = .false. + errmsg = 'at end of line looking for floating point number' + ierrpos=-1 + return + endif + jrec=jrec+1 + i1=istrt(jrec) + i2=i1+inumb(jrec)-1 + ie2=i2 +c... sign + isign=1 + if (ia(i1:i1).eq.xchar(12))isign=-1 + if (ia(i1:i1).eq.xchar(12).or.ia(i1:i1).eq.xchar(11)) i1=i1+1 +c... exponent + do ie=i1,i2 + if (ia(ie:ie).eq.xchar(14) .or. ia(ie:ie).eq.xchar(15)) goto 20 + enddo + iexp=0 + go to 50 + 20 i2=ie-1 + iexp=1 + ie1=ie+1 + if (ia(ie1:ie1).eq.xchar(12))iexp=-1 + if (ia(ie1:ie1).eq.xchar(12).or.ia(ie1:ie1).eq.xchar(11)) + * ie1=ie1+1 + ibuff=0 + do i=ie1,ie2 + do j=1,10 + if (ia(i:i).eq.xchar(j)) go to 41 + enddo + goto 100 + 41 ibuff=ibuff*10+j-1 + enddo + iexp=iexp*ibuff +c.... the number itself + 50 orep=.false. + do i=i1,i2 + if(ia(i:i).ne.xchar(13)) then + do j=1,10 + if (ia(i:i).eq.xchar(j)) go to 70 + enddo + goto 100 + 70 buf=buf*ten+ dfloat(j-1) + else + if(orep)go to 100 + iexp=iexp+i-i2 + orep=.true. + endif + enddo + buf = buf * dfloat(isign) * ten**iexp + inp_f = .true. + return +c + 100 inp_f = .false. + jrec = jrec-1 ! Position to re-read the field + ierrpos = i + errmsg = 'illegal character reading floating point number' +c + end +logical function inp_i(jbuf) + implicit none + + character*1 xchar(12) + integer jbuf + data xchar /'0','1','2','3','4','5','6','7','8','9' + 1 ,'+','-'/ +c +c subroutine for reading integers from the array ia, +c starting at ia(istrt(jrec)) and going on for inumb(jrec)) +c elements. plus signs are ignored, the answer is accumulated +c in jbuf +c + ierrpos = -1 + errmsg = ' ' + jbuf = 0 + if(jrec.ge.jump) then + inp_i = .false. + ierrpos = -1 + errmsg = 'at end of line looking for integer' + return + endif + jrec = jrec + 1 + n = inumb(jrec) + ifact = 1 + ist=istrt(jrec) + nstrt = ist + n - 1 + do i = 1,n + xtemp = ia(nstrt:nstrt) + do j=1,12 + if(xchar(j).eq.xtemp)go to 130 + enddo + goto 120 +c + 130 if(j.ge.11) then + if(nstrt.ne.ist)go to 120 + if(j.ge.12)jbuf=-jbuf + go to 160 + endif + jbuf=jbuf+(j-1)*ifact + ifact = ifact * 10 + nstrt=nstrt-1 + enddo + 160 continue + inp_i = .true. + return +c + 120 ierrpos = nstrt + errmsg = 'illegal character when reading integer' + inp_i = .false. + jrec = jrec-1 + return +c + end +logical function inp_compare(ocase, a, b) + implicit none + logical ocase + character*(*) a, b + integer la, lb, i + character*1 atest, btest + integer inp_strlen + external inp_strlen +c + inp_compare = .false. + la = inp_strlen(a) + lb = inp_strlen(b) + if (la .gt. lb) then + return + else if (ocase) then + inp_compare = a .eq. b(1:la) + return + else + do i = 1, la + atest = a(i:i) + btest = b(i:i) + call inp_lcase(atest) + call inp_lcase(btest) + if (atest.ne.btest) return + enddo + inp_compare = .true. + return + endif +c + end +logical function inp_match(nrec, ocase, test, array, ind) + implicit none + integer nrec, ind + logical ocase, inp_compare + character*(*) test, array(*) + integer i, l, inp_strlen + external inp_compare, inp_strlen +c + l = inp_strlen(test) + inp_match = .false. + ind = -1 +c + do i=1,nrec + if (inp_compare(ocase, test(1:l), array(i))) then + if (inp_match) then + inp_match = .false. ! Ambiguity + ind = 0 + return + else + inp_match = .true. ! First match + ind = i + endif + endif + enddo +c + end +subroutine inp_prev_field() + implicit none + +c + call inp_set_field(max(0,inp_cur_field()-1)) +c + end +integer function inp_strlen(a) + implicit none + + character*(*) a + integer i + integer len + logical ois_ws + intrinsic len + ois_ws(xtest) = (xtest.eq.xblnk .or. xtest.eq.xtab) +c + do i = len(a),1,-1 + if (.not. ois_ws(a(i:i))) goto 10 + enddo +c + 10 inp_strlen = i +c + end +subroutine inp_lcase(string) + implicit none + character*(*) string + intrinsic ichar, len + integer i, length, uca, ucz, lca, shift, test +c + uca = ichar('A') ! MUST be uppercase A + ucz = ichar('Z') ! MUST be uppercase Z + lca = ichar('a') ! MUST be lowercase a + shift = lca - uca + if (shift .eq. 0) + $ call errquit('inp_lcase: check case of program source', 0) +c + length = len(string) + do i = 1, length + test = ichar(string(i:i)) + if (test.ge.uca .and. test.le.ucz) + $ string(i:i) = char(test+shift) + enddo +c + end +logical function inp_search(ocase, z) + implicit none + character*(*) z + logical ocase + character*256 tmp + integer length + integer inp_strlen + logical inp_read, inp_a, inp_compare + external inp_read, inp_a, inp_compare, inp_strlen +c + length = inp_strlen(z) +c + 10 if (inp_read()) then + if (inp_a(tmp)) then + if (inp_compare(ocase, z(1:length), tmp)) then + call inp_prev_field() + inp_search = .true. + return + endif + endif + goto 10 + endif +c + inp_search = .false. +c + end \ No newline at end of file diff --git a/src/inp/inp.h b/src/inp/inp.h index d86cb78..ec1f638 100644 --- a/src/inp/inp.h +++ b/src/inp/inp.h @@ -1,18 +1,19 @@ #ifndef _INP_H #define _INP_H +#include + + bool inp_i(int); + bool inp_f(double); + bool inp_a(char *); + bool inp_read(); + bool inp_line(char *); + bool inp_match(int, bool, char *, char *[], int); + bool inp_search(bool, char*); + bool inp_compare(bool, char *, char *); + bool inp_eof(); + bool inp_a_trunc(char *); - int inp_i(); - float inp_f(); - char *inp_a(); int inp_n_field(); int inp_cur_field(); - int inp_match(); - int inp_read(); - int inp_search(); - int inp_compare(); - int inp_strlen(); - int inp_line(); - int inp_eof(); - int inp_a_trunc(); #endif \ No newline at end of file diff --git a/src/inp/inpP.h b/src/inp/inpP.h index 0ddece8..f1d8a52 100644 --- a/src/inp/inpP.h +++ b/src/inp/inpP.h @@ -7,7 +7,7 @@ #define MAX_WIDTH 256 // Maximum no. of characters in an input line #define MAX_FIELD MAX_WIDTH/2 + 1 // Maximum no. of fields in an input line - char ja[256]; // Input buffers ... MUST match max_width + char ja[256], ia[256]; // Input buffers ... MUST match max_width char tmp[256]; // Same size work space char errmsg[80]; // Error message char xcomm; // Comment character diff --git a/src/input/GNUmakefile b/src/input/GNUmakefile new file mode 100644 index 0000000..a67d17a --- /dev/null +++ b/src/input/GNUmakefile @@ -0,0 +1,8 @@ + + OBJ = input_parse.o input_mem_size.o memory_input.o input_set.o \ + input_start_opt.o input_title.o + LIBRARY = libinput.a + +include ../config/makefile.h +include ../config/makelib.h + diff --git a/src/input/design b/src/input/design new file mode 100644 index 0000000..ac7c3f7 --- /dev/null +++ b/src/input/design @@ -0,0 +1,18 @@ + + + separate routine input_mem_size() scans input for memory directive + + separate routine input_rtdb_name() scans input to infer the rtdb name + + top level recognizes simple directives and module names only + + inp_read() + while (input available) + + read name + + match name against known directives and call appropriate + modules to handle the input + + + diff --git a/src/input/input.format b/src/input/input.format new file mode 100644 index 0000000..f037a2b --- /dev/null +++ b/src/input/input.format @@ -0,0 +1,179 @@ + +1) All input is free format and is lower cased on input except for + file names and titles. + +2) Directive structure + +3) Most directives can appear in any order + +4) Sensible defaults + full error checking +--------------------- + + +--------------------------------------------------------------------- +Directives +---------- + + Syntax for definition of the directives + + () used to group entries (not actually present in the input) + || separate exclusive formats + [] enclose optional entries with a default value + <> enclose a type and a name of a value to be specified + A string is just a sequence of characters, enclosed in + quotes if there is white space + \ is used to concatenate lines + + The order of keyed optional entries should not matter + unless noted otherwise. + +--------------------------------------------------------------------- +The input must commence with either a START or a RESTART directive +which have the same syntax + +(RESTART || START) \ + [[PREFIX] = (-'.db' || 'calc')] \ + [DATABASE = .db] + + These directives determine if this is a restart or startup calculation + and provide definition of and + In a startup calculation any existing data base is destroyed. + + By default all filenames will be created by appending to a common + file prefix, which could include a path adjustment. This defaults + to either the data base name, stripped of a trailing '.db', or + failing that 'calc'. + + The data base path can be specified, or defaulted using the + file prefix. + + E.g. + + start + + Startup using all defaults ( = 'calc') and + = 'calc.db' + + restart water + + Restart calculation with = 'water' and + = 'water.db' + + restart prefix water + + Same as previous example + + restart /tmp/rjh/ch2 database /tmp/rjh/ch2small.db + + Restart calculation with = '/tmp/rjh/ch2' and + = '/tmp/rjh/ch2/ch2small.db' + + start database /disk2/mgo.db + + Startup calculation with = '/disk2/mgo' and + = '/disk2/mgo.db' + + start database /disk2/mgo_dumpfile + + Startup calculation with = 'calc' and + = '/disk2/mgo_dumpfile' + +--------------------------------------------------------------------- + +TITLE + + Enters the string into the data base entry 'title' + +--------------------------------------------------------------------- + +GEOMETRY [ = 'geometry'] [[UNITS] = 'au'] + + read until encounter END + + +END + + Enters atomic cartesian coordinates in either atomic units (units + = 'au') or angstroms (units = 'angstroms'). + + Geometries may be optionally named, however, the default name of + 'geometry' must usually be present for a calculation to proceed. + + The atomic tag serves to match against tags provided for basis + function centers. Also, the first 1 or 2 characters of the atomic + tag may interpreted to identify the element. + + e.g. + + geometry 'water at 90 degrees' angs + o 8 0.0 0.0 0.0 + h 1 1.0 0.0 0.0 + h 1 0.0 1.0 0.0 + end +--------------------------------------------------------------------- + +BASIS [ = 'mo basis set'>] LIBRARY + Read until END encountered + + LIBRARY + + or + + + read until next or END encountered + + +END + + If the basis directive is not provided in a startup calculation + then a default of 3-21g is adopted. + + Basis sets may also be named, with the default name of 'mo basis set' + being that required by modules that compute MOs. + + Many standard basis sets are available in a library. These may be + used for the whole molecule or just for individual atoms. + Basis functions defined within the basis set directive add to those + adopted on the directive line. + + e.g. + + basis library ccpvdz + + Just use the standard cpvdz basis set + + basis + h s + 0.01 100. + 0.8 7. + h p + 1.0 1.0 + end + + Defines a rather stupid basis set for the atomic center with tag h. + + basis library 3-21g + o d + 1.0 0.001 + si library "somebody's standard diffuse polarization functions" + end + + Adopts a 3-1g basis set which is augmented with a d function + on centers with tag o and a standard set of functions on centers + with tag si. + +--------------------------------------------------------------------- + +RHF + + + + +--------------------------------------------------------------------- + +The CALCULATION directive will eventually control what calculations +are performed and in what order, with high level control of input +and output to each module. Right now, since there is only RHF +energy, it does very little ... in fact it is not even implemented. + +--------------------------------------------------------------------- diff --git a/src/util/itri.h b/src/util/itri.h new file mode 100644 index 0000000..e4bf5e6 --- /dev/null +++ b/src/util/itri.h @@ -0,0 +1,13 @@ +#ifndef _ITRI_H +#define _ITRI_H + +/* +c simple statement function for evaluating index into lower +c triangular packed array ... indices do not have to be ordered +c +c This file must be include immediately before the first executable +c statement +*/ + itri[i][j] = (max(i,j)*(max(i,j)-3))/2 + i + j; + +#endif // _ITRI_H \ No newline at end of file diff --git a/src/util/output.c b/src/util/output.c new file mode 100644 index 0000000..cb0b747 --- /dev/null +++ b/src/util/output.c @@ -0,0 +1,290 @@ +#include + +void output(double *z, int rowlow, int rowhi, int collow, int colhi, + int rowdim, int coldim, int nctl) { + +/*...................................................................... +c output prints a real*8 matrix in formatted form with numbered rows +c and columns. the input is as follows; +c matrix(*,*).........matrix to be output +c rowlow..............row number at which output is to begin +c rowhi...............row number at which output is to end +c collow..............column number at which output is to begin +c colhi...............column number at which output is to end +c rowdim..............row dimension of matrix(*,*) +c coldim..............column dimension of matrix(*,*) +c nctl................carriage control flag; 1 for single space +c 2 for double space +c 3 for triple space +c the parameters that follow matrix are all of type integer*4. the +c program is set up to handle 5 columns/page with a 1p5d24.15 format for +c the columns. if a different number of columns is required, change +c formats 1000 and 2000, and initialize kcol with the new number of +c columns. +c author; nelson h.f. beebe, quantum theory project, university of +c florida, gainesville +c....................................................................... +C$Id$*/ + int begin, kcol, nctl, i, j, last, k; + double zero = 0.0; + char asa[3][8] = {" ", "00000000", "--------"}; + char ctl, blank = ' '; + + kcol = 8; + if (rowhi < rowlow || colhi < collow) { + printf(" zero matrix\n"); + return; + } + + last = (colhi < collow + kcol - 1) ? colhi : collow + kcol - 1; + for (begin = collow; begin <= colhi; begin += kcol) { + for (i = begin; i <= last; i++) { + printf("%d ", i); + } + printf("\n"); + for (k = rowlow; k <= rowhi; k++) { + for (i = begin; i <= last; i++) { + if (z[k * rowdim + i] != zero) { + printf("%d ", k); + for (i = begin; i <= last; i++) { + printf("%f ", z[k * rowdim + i]); + } + printf("\n"); + break; + } + } + } + last = (last + kcol < colhi) ? last + kcol : colhi; + } +} + + subroutine zoutput (z,rowlow,rowhi,collow,colhi,rowdim,coldim, + $ nctl) +c....................................................................... +c output prints a complex*16 matrix in formatted form with numbered rows +c and columns. the input is as follows; +c matrix(*,*).........matrix to be output +c rowlow..............row number at which output is to begin +c rowhi...............row number at which output is to end +c collow..............column number at which output is to begin +c colhi...............column number at which output is to end +c rowdim..............row dimension of matrix(*,*) +c coldim..............column dimension of matrix(*,*) +c nctl................carriage control flag; 1 for single space +c 2 for double space +c 3 for triple space +c the parameters that follow matrix are all of type integer*4. the +c program is set up to handle 5 columns/page with a 1p5d24.15 format for +c the columns. if a different number of columns is required, change +c formats 1000 and 2000, and initialize kcol with the new number of +c columns. +c author; nelson h.f. beebe, quantum theory project, university of +c florida, gainesville +c....................................................................... +C$Id$ + implicit none + integer rowlow,rowhi,collow,colhi,rowdim,coldim,begin,kcol + integer nctl, i, j, last, k + double complex z(rowdim,coldim), zero + character*8 asa(3), ctl, blank +* character*8 column +* data column/'column' / + data asa/' ','00000000' , + 1 '--------' /,blank/' '/ + data kcol/8/ + data zero/0.d00/ + do 11 i=rowlow,rowhi + do 10 j=collow,colhi + if (z(i,j).ne.zero) go to 15 + 10 continue + 11 continue + write (6,3000) + 3000 format (/' zero matrix'/) + go to 3 + 15 continue + ctl = blank + if ((nctl.le.3).and.(nctl.gt.0)) ctl = asa(nctl) + if (rowhi.lt.rowlow) go to 3 + if (colhi.lt.collow) go to 3 + last = min(colhi,collow+kcol-1) + do 2 begin = collow,colhi,kcol +* write (6,1000) (column,i,i = begin,last) + write (6,1000) (i,i = begin,last) + do 1 k = rowlow,rowhi + do 4 i=begin,last + if (z(k,i).ne.zero) go to 5 + 4 continue + go to 1 + 5 write (6,2000) ctl,k,(z(k,i), i = begin,last) + 1 continue + last = min(last+kcol,colhi) + 2 continue + 3 return +* kcol = 4 +* 1000 format (/1h ,16x,3(a6,i3,2x),(a6,i3)) +* 2000 format (a1,3hrow,i4,2x,4f17.11) +* kcol = 8 +* +* if U like having rows and columns labelled with row and col +* use these +* +* 1000 format (/1h ,11x,7(a3,i3,3x),(a3,i3)) +* 2000 format (a1,'row',i4,1x,8f9.4) +c + 1000 format (/1h ,8x,7(' ',i3,3x),(' ',i3)) + 2000 format (a1,i4,1x,1p,16d9.2) + end +c + subroutine doutput (z,rowlow,rowhi,collow,colhi,rowdim,coldim, + $ nctl) +c....................................................................... +c output prints a real*8 matrix in formatted form with numbered rows +c and columns. the input is as follows; +c matrix(*,*).........matrix to be output +c rowlow..............row number at which output is to begin +c rowhi...............row number at which output is to end +c collow..............column number at which output is to begin +c colhi...............column number at which output is to end +c rowdim..............row dimension of matrix(*,*) +c coldim..............column dimension of matrix(*,*) +c nctl................carriage control flag; 1 for single space +c 2 for double space +c 3 for triple space +c the parameters that follow matrix are all of type integer*4. the +c program is set up to handle 5 columns/page with a 1p5d24.15 format for +c the columns. if a different number of columns is required, change +c formats 1000 and 2000, and initialize kcol with the new number of +c columns. +c author; nelson h.f. beebe, quantum theory project, university of +c florida, gainesville +c....................................................................... +C$Id$ + implicit none + integer rowlow,rowhi,collow,colhi,rowdim,coldim,begin,kcol + integer nctl, i, j, last, k + double precision z(rowdim,coldim), zero + character*8 asa(3), ctl, blank +* character*8 column +* data column/'column' / + data asa/' ','00000000' , + 1 '--------' /,blank/' '/ + data kcol/8/ + data zero/0.d00/ + do 11 i=rowlow,rowhi + do 10 j=collow,colhi + if (z(i,j).ne.zero) go to 15 + 10 continue + 11 continue + write (6,3000) + 3000 format (/' zero matrix'/) + go to 3 + 15 continue + ctl = blank + if ((nctl.le.3).and.(nctl.gt.0)) ctl = asa(nctl) + if (rowhi.lt.rowlow) go to 3 + if (colhi.lt.collow) go to 3 + last = min(colhi,collow+kcol-1) + do 2 begin = collow,colhi,kcol +* write (6,1000) (column,i,i = begin,last) + write (6,1000) (i,i = begin,last) + do 1 k = rowlow,rowhi + do 4 i=begin,last + if (z(k,i).ne.zero) go to 5 + 4 continue + go to 1 + 5 write (6,2000) ctl,k,(z(k,i), i = begin,last) + 1 continue + last = min(last+kcol,colhi) + 2 continue + 3 return +* kcol = 4 +* 1000 format (/1h ,16x,3(a6,i3,2x),(a6,i3)) +* 2000 format (a1,3hrow,i4,2x,4f17.11) +* kcol = 8 +* +* if U like having rows and columns labelled with row and col +* use these +* +* 1000 format (/1h ,11x,7(a3,i3,3x),(a3,i3)) +* 2000 format (a1,'row',i4,1x,8f9.4) +c + 1000 format (/1h ,8x,7(' ',i3,3x),(' ',i3)) + 2000 format (a1,i4,1x,1p,8d9.2) + end +c + subroutine ioutput (z,rowlow,rowhi,collow,colhi,rowdim,coldim, + $ nctl) +c....................................................................... +c output prints a real*8 matrix in formatted form with numbered rows +c and columns. the input is as follows; +c matrix(*,*).........matrix to be output +c rowlow..............row number at which output is to begin +c rowhi...............row number at which output is to end +c collow..............column number at which output is to begin +c colhi...............column number at which output is to end +c rowdim..............row dimension of matrix(*,*) +c coldim..............column dimension of matrix(*,*) +c nctl................carriage control flag; 1 for single space +c 2 for double space +c 3 for triple space +c the parameters that follow matrix are all of type integer*4. the +c program is set up to handle 5 columns/page with a 1p5d24.15 format for +c the columns. if a different number of columns is required, change +c formats 1000 and 2000, and initialize kcol with the new number of +c columns. +c author; nelson h.f. beebe, quantum theory project, university of +c florida, gainesville +c....................................................................... +C$Id$ + implicit none + integer rowlow,rowhi,collow,colhi,rowdim,coldim,begin,kcol + integer nctl, i, j, last, k + integer z(rowdim,coldim), zero + character*8 asa(3), ctl, blank +* character*8 column +* data column/'column' / + data asa/' ','00000000' , + 1 '--------' /,blank/' '/ + data kcol/8/ + data zero/0/ + do 11 i=rowlow,rowhi + do 10 j=collow,colhi + if (z(i,j).ne.zero) go to 15 + 10 continue + 11 continue + write (6,3000) + 3000 format (/' zero matrix'/) + go to 3 + 15 continue + ctl = blank + if ((nctl.le.3).and.(nctl.gt.0)) ctl = asa(nctl) + if (rowhi.lt.rowlow) go to 3 + if (colhi.lt.collow) go to 3 + last = min(colhi,collow+kcol-1) + do 2 begin = collow,colhi,kcol +* write (6,1000) (column,i,i = begin,last) + write (6,1000) (i,i = begin,last) + do 1 k = rowlow,rowhi + do 4 i=begin,last + if (z(k,i).ne.zero) go to 5 + 4 continue + go to 1 + 5 write (6,2000) ctl,k,(z(k,i), i = begin,last) + 1 continue + last = min(last+kcol,colhi) + 2 continue + 3 return +* kcol = 4 +* 1000 format (/1h ,16x,3(a6,i3,2x),(a6,i3)) +* 2000 format (a1,3hrow,i4,2x,4f17.11) +* kcol = 8 +* +* if U like having rows and columns labelled with row and col +* use these +* +* 1000 format (/1h ,11x,7(a3,i3,3x),(a3,i3)) +* 2000 format (a1,'row',i4,1x,8f9.4) +c + 1000 format (/1h ,8x,7(' ',i3,3x),(' ',i3)) + 2000 format (a1,i4,1x,8i9) + end \ No newline at end of file diff --git a/src/util/util.h b/src/util/util.h new file mode 100644 index 0000000..50c45a0 --- /dev/null +++ b/src/util/util.h @@ -0,0 +1,6 @@ +#ifndef _UTIL_H +#define _UTIL_H + +double ddot(); + +#endif // _UTIL_H \ No newline at end of file