slasq5.c

最著名最快的分子模拟软件

#include <math.h>
#include "gmx_lapack.h"

void
F77_FUNC(slasq5,SLASQ5)(int *i0, 
	int *n0,
	float *z__, 
	int *pp, 
	float *tau,
	float *dmin__, 
	float *dmin1, 
	float *dmin2, 
	float *dn,
	float *dnm1, 
	float *dnm2,
	int *ieee)
{
    int i__1;
    float d__1, d__2;

    static float d__;
    static int j4, j4p2;
    static float emin, temp;

    --z__;

    if (*n0 - *i0 - 1 <= 0) {
	return;
    }

    j4 = (*i0 << 2) + *pp - 3;
    emin = z__[j4 + 4];
    d__ = z__[j4] - *tau;
    *dmin__ = d__;
    *dmin1 = -z__[j4];

    if (*ieee) {

	if (*pp == 0) {
	    i__1 = 4*(*n0 - 3);
	    for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
		z__[j4 - 2] = d__ + z__[j4 - 1];
		temp = z__[j4 + 1] / z__[j4 - 2];
		d__ = d__ * temp - *tau;
                if(d__<*dmin__)
                  *dmin__ = d__;
		z__[j4] = z__[j4 - 1] * temp;
		d__1 = z__[j4];
                if(d__1<emin)
                  emin = d__1;
	    }
	} else {
	    i__1 = 4*(*n0 - 3);
	    for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
		z__[j4 - 3] = d__ + z__[j4];
		temp = z__[j4 + 2] / z__[j4 - 3];
		d__ = d__ * temp - *tau;
                if(d__<*dmin__)
                  *dmin__ = d__;
		z__[j4 - 1] = z__[j4] * temp;
		d__1 = z__[j4 - 1];
                if(d__1<emin)
                  emin = d__1;
	    }
	}

	*dnm2 = d__;
	*dmin2 = *dmin__;
	j4 = 4*(*n0 - 2) - *pp;
	j4p2 = j4 + (*pp << 1) - 1;
	z__[j4 - 2] = *dnm2 + z__[j4p2];
	z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
	*dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau;
        if(*dnm1<*dmin__)
          *dmin__ = *dnm1;

	*dmin1 = *dmin__;
	j4 += 4;
	j4p2 = j4 + (*pp << 1) - 1;
	z__[j4 - 2] = *dnm1 + z__[j4p2];
	z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
	*dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau;
        if(*dn<*dmin__)
          *dmin__ = *dn;

    } else {

	if (*pp == 0) {
	    i__1 = 4*(*n0 - 3);
	    for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
		z__[j4 - 2] = d__ + z__[j4 - 1];
		if (d__ < 0.) {
		    return;
		} else {
		    z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
		    d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]) - *tau;
		}
                if(d__<*dmin__)
                  *dmin__ = d__;
		d__1 = emin, d__2 = z__[j4];
		emin = (d__1<d__2) ? d__1 : d__2;
	    }
	} else {
	    i__1 = 4*(*n0 - 3);
	    for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
		z__[j4 - 3] = d__ + z__[j4];
		if (d__ < 0.) {
		    return;
		} else {
		    z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
		    d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]) - *tau;
		}
                if(d__<*dmin__)
                  *dmin__ = d__;
		d__1 = emin, d__2 = z__[j4 - 1];
		emin = (d__1<d__2) ? d__1 : d__2;
	    }
	}

	*dnm2 = d__;
	*dmin2 = *dmin__;
	j4 = 4*(*n0 - 2) - *pp;
	j4p2 = j4 + (*pp << 1) - 1;
	z__[j4 - 2] = *dnm2 + z__[j4p2];
	if (*dnm2 < 0.) {
	    return;
	} else {
	    z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
	    *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau;
	}
        if(*dnm1<*dmin__)
          *dmin__ = *dnm1;

	*dmin1 = *dmin__;
	j4 += 4;
	j4p2 = j4 + (*pp << 1) - 1;
	z__[j4 - 2] = *dnm1 + z__[j4p2];
	if (*dnm1 < 0.) {
	    return;
	} else {
	    z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
	    *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau;
	}
        if(*dn<*dmin__)
          *dmin__ = *dn;

    }

    z__[j4 + 2] = *dn;
    z__[(*n0 << 2) - *pp] = emin;
    return;

}