continued translation of files (mostly inp.c)

This commit is contained in:
Adam Parler 2024-01-21 00:01:51 -08:00
parent ec4e4388fe
commit d4bd5d13bf
11 changed files with 1159 additions and 22 deletions

View file

@ -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

View file

@ -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

630
src/inp/inp.c Normal file
View file

@ -0,0 +1,630 @@
#include <string.h>
#include <stdbool.h>
#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 <stdbool.h>
#include <string.h>
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

View file

@ -1,18 +1,19 @@
#ifndef _INP_H
#define _INP_H
#include <stdbool.h>
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

View file

@ -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

8
src/input/GNUmakefile Normal file
View file

@ -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

18
src/input/design Normal file
View file

@ -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

179
src/input/input.format Normal file
View file

@ -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] <string file prefix> = (<data base path>-'.db' || 'calc')] \
[DATABASE <string data base path> = <file prefix>.db]
These directives determine if this is a restart or startup calculation
and provide definition of <file prefix> and <data base path>
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 (<file prefix> = 'calc') and
<data base path> = 'calc.db'
restart water
Restart calculation with <file prefix> = 'water' and
<data base path> = 'water.db'
restart prefix water
Same as previous example
restart /tmp/rjh/ch2 database /tmp/rjh/ch2small.db
Restart calculation with <file prefix> = '/tmp/rjh/ch2' and
<data base path> = '/tmp/rjh/ch2/ch2small.db'
start database /disk2/mgo.db
Startup calculation with <file prefix> = '/disk2/mgo' and
<data base path> = '/disk2/mgo.db'
start database /disk2/mgo_dumpfile
Startup calculation with <file prefix> = 'calc' and
<data base path> = '/disk2/mgo_dumpfile'
---------------------------------------------------------------------
TITLE <string title>
Enters the string into the data base entry 'title'
---------------------------------------------------------------------
GEOMETRY [<string name> = 'geometry'] [[UNITS] <string units> = 'au']
read until encounter END
<string atom tag> <double charge> <double x> <double y> <double z>
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 [<string name> = 'mo basis set'>] LIBRARY <string standard set>
Read until END encountered
<string atom tag> LIBRARY <string standard set>
or
<string atom tag> <string angular momentum>
read until next <atom tag> or END encountered
<double contraction coeff> <double exponent>
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.
---------------------------------------------------------------------

13
src/util/itri.h Normal file
View file

@ -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

290
src/util/output.c Normal file
View file

@ -0,0 +1,290 @@
#include <stdio.h>
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

6
src/util/util.h Normal file
View file

@ -0,0 +1,6 @@
#ifndef _UTIL_H
#define _UTIL_H
double ddot();
#endif // _UTIL_H