/* PSolve.f -- translated by f2c (version of 20 August 1993  13:15:44).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    char curpform[5];
} forms_;

#define forms_1 forms_

struct {
    integer n, lda;
} matdim_;

#define matdim_1 matdim_

struct {
    real a[40000], m[200];
} system_;

#define system_1 system_

/* Table of constant values */

static integer c__9 = 9;
static integer c__1 = 1;


/*  This file contains the preconditioner solve routines: */

/*  PSOLVE and PSOLVETRANS call the appropriate solver: */

/*     PSOLVENONE and PSOLVENONETRANS for using no preconditioning. */

/*     PSOLVEJAC and PSOLVEJACTRANS for Jacobi preconditioning. */

/*     Also included are the solvers for QMR which require left and right */
/*     preconditioning: PSOLVEQ and PSOLVETRANSQ */

/*     ================================================================ */
/* Subroutine */ int psolve_(x, b)
real *x, *b;
{
    /* Builtin functions */
    integer s_wsle(), do_lio(), e_wsle();
    /* Subroutine */ int s_stop();

    /* Local variables */
    extern logical lsame_();
    extern /* Subroutine */ int psolvejac_(), psolvenone_();

    /* Fortran I/O blocks */
    static cilist io___1 = { 0, 6, 0, 0, 0 };



/*     .. Array Arguments .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Common Blocks .. */

/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    if (lsame_(forms_1.curpform, "IDENT", 5L, 5L)) {
	psolvenone_(&x[1], &b[1]);
    } else if (lsame_(forms_1.curpform, "JACBI", 5L, 5L)) {
	psolvejac_(&x[1], &b[1]);
    } else {
	s_wsle(&io___1);
	do_lio(&c__9, &c__1, "IN PSOLVE: UNKNOWN PRECONDITIONER", 33L);
	do_lio(&c__9, &c__1, forms_1.curpform, 5L);
	do_lio(&c__9, &c__1, " QUITTING", 9L);
	e_wsle();
	s_stop("", 0L);
    }

    return 0;

} /* psolve_ */

/*     ================================================================ */

/* Subroutine */ int psolvetrans_(x, b)
real *x, *b;
{
    /* Builtin functions */
    integer s_wsle(), do_lio(), e_wsle();
    /* Subroutine */ int s_stop();

    /* Local variables */
    extern /* Subroutine */ int psolvenonetrans_();
    extern logical lsame_();
    extern /* Subroutine */ int psolvejactrans_();

    /* Fortran I/O blocks */
    static cilist io___2 = { 0, 6, 0, 0, 0 };



/*     .. Array Arguments .. */
/*     .. */
/*     .. Common Blocks .. */
/*     .. */
/*     .. Local Scalars .. */

/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    if (lsame_(forms_1.curpform, "IDENT", 5L, 5L)) {
	psolvenonetrans_(&x[1], &b[1]);
    } else if (lsame_(forms_1.curpform, "JACBI", 5L, 5L)) {
	psolvejactrans_(&x[1], &b[1]);
    } else {
	s_wsle(&io___2);
	do_lio(&c__9, &c__1, "IN PSOLVE: UNKNOWN PRECONDITIONER", 33L);
	do_lio(&c__9, &c__1, forms_1.curpform, 5L);
	do_lio(&c__9, &c__1, " QUITTING", 9L);
	e_wsle();
	s_stop("", 0L);
    }

    return 0;

} /* psolvetrans_ */

/*     ================================================================ */

/* Subroutine */ int psolvenone_(x, b)
real *x, *b;
{
    extern /* Subroutine */ int scopy_();


/*     .. Array Arguments .. */
/*     .. */
/*     .. Common Blocks .. */


/*  Purpose */
/*  ======= */

/*  This PSOLVE is for the unpreconditioned version, i.e. just does */
/*  a vector copy ( B to X ) then returns. */

/*  Arguments */
/*  ========= */

/*  B       (input) REAL             array, dimension N. */
/*          On entry, right hand side vector B. */
/*          Unchanged on exit. */

/*  X       (output) REAL             array, dimension N. */
/*          Set to solution on output. */

/*  BLAS:  SCOPY */
/*  ============================================================ */

/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    scopy_(&matdim_1.n, &b[1], &c__1, &x[1], &c__1);

    return 0;

/*     End of PSolveNone */

} /* psolvenone_ */


/*     ===================================================== */
/* Subroutine */ int psolvenonetrans_(x, b)
real *x, *b;
{
    extern /* Subroutine */ int scopy_();


/*     .. Array Arguments .. */
/*     .. */
/*     .. */
/*     .. Common Blocks .. */

/*     .. */

/*  Purpose */
/*  ======= */

/*  This PSOLVE is for the unpreconditioned version, i.e. just does */
/*  a vector copy ( B to X ) then returns. */

/*  Arguments */
/*  ========= */

/*  B       (input) REAL             array, dimension N. */
/*          On entry, right hand side vector B. */
/*          Unchanged on exit. */

/*  X       (output) REAL             array, dimension N. */
/*          Set to solution on output. */

/*  BLAS:  SCOPY */
/*  ============================================================ */

/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    scopy_(&matdim_1.n, &b[1], &c__1, &x[1], &c__1);

    return 0;

/*     End of PSolve */

} /* psolvenonetrans_ */


/*     =========================================================== */
/* Subroutine */ int psolvejac_(x, b)
real *x, *b;
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i;


/*     .. Array Delcarations .. */
/*     .. */
/*     .. Common Blocks .. */




/*  Purpose */
/*  ======= */

/*  PSOLVE solves the linear system Mx = b where matrix M has */
/*  is diagonal. */

/*  Arguments */
/*  ========= */

/*  B       (input) REAL             array, dimension N. */
/*          On entry, right hand side vector B. */
/*          Unchanged on exit. */

/*  X       (output) REAL             array, dimension N. */
/*          Set to solution on output. */
/*  ============================================================ */

/*     .. Local Scalars .. */


/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    i__1 = matdim_1.n;
    for (i = 1; i <= i__1; ++i) {
	x[i] = b[i] / system_1.m[i - 1];
/* L10: */
    }

    return 0;

/*     End of PSolveJac */

} /* psolvejac_ */


/*     ========================================================= */
/* Subroutine */ int psolvejactrans_(x, b)
real *x, *b;
{
    extern /* Subroutine */ int psolvejac_();


/*     .. Array Delcarations .. */

/*     .. */

/*  Purpose */
/*  ======= */

/*  PSOLVETRANS solves the linear system Mx = b where matrix M has */
/*  is diagonal. Since this is the same as the non-transpose version, */
/*  this routine is actual just a mask to PSOLVEJAC. */

/*  Arguments */
/*  ========= */

/*  B       (input) REAL             array, dimension N. */
/*          On entry, right hand side vector B. */
/*          Unchanged on exit. */

/*  X       (output) REAL             array, dimension N. */
/*          Set to solution on output. */
/*  ============================================================ */

/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    psolvejac_(&x[1], &b[1]);

    return 0;

/*     End of PSolveJacTrans */

} /* psolvejactrans_ */

/*     ================================================================ */
/*     Following are the solvers for QMR, allowing left and right */
/*     preconditioning. */
/*     ================================================================ */
/* Subroutine */ int psolveq_(x, b, which, which_len)
real *x, *b;
char *which;
ftnlen which_len;
{
    /* Builtin functions */
    integer s_wsle(), do_lio(), e_wsle();
    /* Subroutine */ int s_stop();

    /* Local variables */
    extern logical lsame_();
    extern /* Subroutine */ int psolvejac_(), psolvenone_();

    /* Fortran I/O blocks */
    static cilist io___4 = { 0, 6, 0, 0, 0 };



/*     .. Argument Declarations .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */

/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    if (lsame_(forms_1.curpform, "IDENT", 5L, 5L)) {
	psolvenone_(&x[1], &b[1]);
    } else if (lsame_(forms_1.curpform, "JACBI", 5L, 5L)) {
	if (lsame_(which, "LEFT", 4L, 4L)) {
	    psolvejac_(&x[1], &b[1]);
	} else {
	    psolvenone_(&x[1], &b[1]);
	}
    } else {
	s_wsle(&io___4);
	do_lio(&c__9, &c__1, "IN PSOLVEQ: UNKNOWN PRECONDITIONER", 34L);
	do_lio(&c__9, &c__1, forms_1.curpform, 5L);
	do_lio(&c__9, &c__1, " QUITTING", 9L);
	e_wsle();
	s_stop("", 0L);
    }

    return 0;

} /* psolveq_ */

/*     ================================================================ */

/* Subroutine */ int psolvetransq_(x, b, which, which_len)
real *x, *b;
char *which;
ftnlen which_len;
{
    /* Builtin functions */
    integer s_wsle(), do_lio(), e_wsle();
    /* Subroutine */ int s_stop();

    /* Local variables */
    extern logical lsame_();
    extern /* Subroutine */ int psolvejac_(), psolvenone_();

    /* Fortran I/O blocks */
    static cilist io___5 = { 0, 6, 0, 0, 0 };



/*     .. Argument Declarations .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Common Blocks .. */

/*     .. Executable Statements .. */

    /* Parameter adjustments */
    --b;
    --x;

    /* Function Body */
    if (lsame_(forms_1.curpform, "IDENT", 5L, 5L)) {
	psolvenone_(&x[1], &b[1]);
    } else if (lsame_(forms_1.curpform, "JACBI", 5L, 5L)) {
	if (lsame_(which, "LEFT", 4L, 4L)) {
	    psolvejac_(&x[1], &b[1]);
	} else {
	    psolvenone_(&x[1], &b[1]);
	}
    } else {
	s_wsle(&io___5);
	do_lio(&c__9, &c__1, "IN PSOLVEQ: UNKNOWN PRECONDITIONER", 34L);
	do_lio(&c__9, &c__1, forms_1.curpform, 5L);
	do_lio(&c__9, &c__1, " QUITTING", 9L);
	e_wsle();
	s_stop("", 0L);
    }

    return 0;

} /* psolvetransq_ */

