/*****************************
 *  FUNCTION:  file2collvec  *
 *****************************/
void file2collvec (V, filename, xm, ym, num, loc, glob)
struct collvec *V;
char filename[80];
int xm, ym, num, loc, glob;
{
int i, j, k, l;
FILE *fp;

fp = fopen(filename, "r");

/*
if ((xm<ym && num==0) || (xm>=ym && num==1))   v2cvA 
    {
    for (i = 0; i < loc; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &V->F[0].el[i].el[j]);

    for (i = 0; i <= glob-2; i++)
        for (l = 0; l < 2; l++)
	    for (j = 0; j < loc; j++)
	        for (k = 0; k < 2; k++)
		    fscanf(fp, "%lf", &V->M[i].el[j].el[k].el[l]);

    for (i = 0; i < loc; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &V->L[0].el[i].el[j]);
    }
else                                           v2cvB 
    {
*/
    for (l = 0; l < loc; l++)
        for (i = 0; i < 2; i++)
	    {
	    fscanf(fp, "%lf", &V->F[0].el[l].el[i]);
	    for (j = 0; j <= glob-2; j++)
	        for (k = 0; k < 2; k++)
		    fscanf(fp, "%lf", &V->M[j].el[l].el[i].el[k]);
	    fscanf(fp, "%lf", &V->L[0].el[l].el[i]);
	    }
/*
    }
*/

fclose(fp);
}  /* end of function file2collvec */


/****************************************************************
 *  FUNCTION:  load_SW_el -- loads the SW element of a collmat  *
 *  M = collmat                                                 *
 *  v = vector from which the collmat is loaded                 *
 *  f = 0, use xl_yg numbering                                  *
 *  f = 1, use yl_xg numbering                                  *
 ****************************************************************/
void load_SW_el (M, fp, f)
struct collmat *M;
FILE *fp;
int f;
{
int i, j, k;

if (f == 0)  /* xl_yg */
    {
    /*  SW equation  */
    fscanf(fp, "%lf", &M->AF[0].a[0].el[0][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->BF[0].a[0].el[0].el[0][i]);
    fscanf(fp, "%lf", &M->AF[0].a[0].el[0][1]);
    fscanf(fp, "%lf", &M->AF[0].b[0].el[0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].a[0].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->BF[0].b[0]      .el[0][i]);
        }

    /*  NW equation  */
    fscanf(fp, "%lf", &M->CF[0].a[0].el[0][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[0].a[0].A.el[0][i]);
    fscanf(fp, "%lf", &M->CF[0].a[0].el[0][1]);
    fscanf(fp, "%lf", &M->CF[0].b[0].el[0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].a[0].B    .el[0][i]);
	fscanf(fp, "%lf", &M->A[0].b[0].el[0].el[0][i]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M->AF[0].a[0].el[1][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->BF[0].a[0].el[0].el[1][i]);
    fscanf(fp, "%lf", &M->AF[0].a[0].el[1][1]);
    fscanf(fp, "%lf", &M->AF[0].b[0].el[1]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].a[0].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->BF[0].b[0]      .el[1][i]);
        }

    /*  NE equation  */
    fscanf(fp, "%lf", &M->CF[0].a[0].el[1][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[0].a[0].C.el[0][i]);
    fscanf(fp, "%lf", &M->CF[0].a[0].el[1][1]);
    fscanf(fp, "%lf", &M->CF[0].b[0].el[1]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].a[0].D    .el[0][i]);
	fscanf(fp, "%lf", &M->A[0].b[0].el[1].el[0][i]);
        }
    }
else  /* f == 1, yl_xg */
    {
    /*  W equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->AF[0].a[0].el[i][j]);
	fscanf(fp, "%lf", &M->AF[0].b[0].el[i]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M->BF[0].a[0].el[j].el[i][k]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->BF[0].b[0].el[i][j]);
        }

    /*  SE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->CF[0].a[0].el[0][j]);
    fscanf(fp, "%lf", &M->CF[0].b[0].el[0]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[0].A.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[0].B.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].b[0].el[0].el[0][j]);

    /*  NE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->CF[0].a[0].el[1][j]);
    fscanf(fp, "%lf", &M->CF[0].b[0].el[1]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[0].C.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[0].D.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].b[0].el[1].el[0][j]);
    }
}  /* end of function load_SW_el */


/**************************************************************
 *  FUNCTION:  load_S_el -- loads the S element of a collmat  *
 *  M = collmat                                               *
 *  v = vector from which the collmat is loaded               *
 *  f = 0, use xl_yg numbering                                *
 *  f = 1, use yl_xg numbering                                *
 *  m = index for which S element                             *
 **************************************************************/
void load_S_el (M, fp, f, m)
struct collmat *M;
FILE *fp;
int f, m;
{
int i, j, mm1, mp1;

if (f == 0)  /* xl_yg */
    {
    mm1 = m - 1;

    /*  SW equation  */
    fscanf(fp, "%lf", &M->AF[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->AF[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->BF[0].a[m]  .el[0].el[0][i]);
        }
    fscanf(fp, "%lf", &M->AF[0].a[m].el[0][1]);
    fscanf(fp, "%lf", &M->AF[0].b[m].el[0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].a[m].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->BF[0].b[m]      .el[0][i]);
        }

    /*  NW equation  */
    fscanf(fp, "%lf", &M->CF[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->CF[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].c[mm1].el[0].el[0][i]);
	fscanf(fp, "%lf", &M->A[0].a[m].A.el[0][i]);
        }
    fscanf(fp, "%lf", &M->CF[0].a[m].el[0][1]);
    fscanf(fp, "%lf", &M->CF[0].b[m].el[0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].a[m].B    .el[0][i]);
	fscanf(fp, "%lf", &M->A[0].b[m].el[0].el[0][i]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M->AF[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->AF[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->BF[0].a[m]  .el[0].el[1][i]);
        }
    fscanf(fp, "%lf", &M->AF[0].a[m].el[1][1]);
    fscanf(fp, "%lf", &M->AF[0].b[m].el[1]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].a[m].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->BF[0].b[m]      .el[1][i]);
        }

    /*  NE equation  */
    fscanf(fp, "%lf", &M->CF[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->CF[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].c[mm1].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->A[0].a[m].C.el[0][i]);
        }
    fscanf(fp, "%lf", &M->CF[0].a[m].el[1][1]);
    fscanf(fp, "%lf", &M->CF[0].b[m].el[1]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].a[m].D    .el[0][i]);
	fscanf(fp, "%lf", &M->A[0].b[m].el[1].el[0][i]);
        }
    }
else  /* f == 1, yl_xg */
    {
    mp1 =m+1;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[0].A.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[0].B.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[0].el[0].el[1][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->B[m].a[0].el[i].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].b[0].el[0][j]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[0].C.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[0].D.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[0].el[1].el[1][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->B[m].a[0].el[i].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].b[0].el[1][j]);

    /*  SE equation  */
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->C[m].a[0].el[i].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].b[0].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[0].A.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[0].B.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].b[0].el[0].el[0][j]);

    /*  NE equation  */
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->C[m].a[0].el[i].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].b[0].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[0].C.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[0].D.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].b[0].el[1].el[0][j]);
    }
}  /* end of function load_S_el */


/****************************************************************
 *  FUNCTION:  load_SE_el -- loads the SE element of a collmat  *
 *  M = collmat                                                 *
 *  v = vector from which the collmat is loaded                 *
 *  f = 0, use xl_yg numbering                                  *
 *  f = 1, use yl_xg numbering                                  *
 *  m = index for which S element                               *
 ****************************************************************/
void load_SE_el (M, fp, f, m)
struct collmat *M;
FILE *fp;
int f, m;
{
int i, j, k, mm1;

if (f == 0)  /* xl_yg */
    {
    mm1 = m - 1;

    /*  SW equation  */
    fscanf(fp, "%lf", &M->AF[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->AF[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->BF[0].a[m]  .el[0].el[0][i]);
        }
    fscanf(fp, "%lf", &M->AF[0].a[m].el[0][1]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->BF[0].a[m].el[1].el[0][i]);
    
    /*  NW equation  */
    fscanf(fp, "%lf", &M->CF[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->CF[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].c[mm1].el[0].el[0][i]);
	fscanf(fp, "%lf", &M->A[0].a[m].A.el[0][i]);
        }
    fscanf(fp, "%lf", &M->CF[0].a[m].el[0][1]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[0].a[m].B    .el[0][i]);

    /*  SE equation  */
    fscanf(fp, "%lf", &M->AF[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->AF[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->BF[0].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->BF[0].a[m]  .el[0].el[1][i]);
        }
    fscanf(fp, "%lf", &M->AF[0].a[m].el[1][1]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->BF[0].a[m].el[1].el[1][i]);

    /*  NE equation  */
    fscanf(fp, "%lf", &M->CF[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->CF[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[0].c[mm1].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->A[0].a[m].C.el[0][i]);
        }
    fscanf(fp, "%lf", &M->CF[0].a[m].el[1][1]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[0].a[m].D.el[0][i]);
    }
else  /* f == 1, yl_xg */
    {
    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[0].A.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[0].B.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[0].el[0].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->BL[0].a[0].el[0][j]);
    fscanf(fp, "%lf", &M->BL[0].b[0].el[0]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[0].C.el[1][j]);
    for (j = 0; j < 2; j++)
	fscanf(fp, "%lf", &M->A[m].a[0].D.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[0].el[1].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->BL[0].a[0].el[1][j]);
    fscanf(fp, "%lf", &M->BL[0].b[0].el[1]);

    /*  E equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M->CL[0].a[0].el[j].el[i][k]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->CL[0].b[0].el[i][j]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->AL[0].a[0].el[i][j]);
	fscanf(fp, "%lf", &M->AL[0].b[0].el[i]);
        }
    }
}  /* end of function load_SE_el */


/**************************************************************
 *  FUNCTION:  load_W_el -- loads the W element of a collmat  *
 *  M = collmat                                               *
 *  v = vector from which the collmat is loaded               *
 *  f = 0, use xl_yg numbering                                *
 *  f = 1, use yl_xg numbering                                *
 *  n = index for which W element                             *
 **************************************************************/
void load_W_el (M, fp, f, n)
struct collmat *M;
FILE *fp;
int f, n;
{
int i, j, k, nm1, np1;

if (f == 0)  /* xl_yg */
    {
    np1 = n + 1;
    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[0].A.el[1][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->B[n].a[0].el[0].el[0][i]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[0].B    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[0].el[0].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].a[0].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->B[n].b[0]      .el[0][i]);
        }

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->C[n].a[0].el[0].el[0][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[np1].a[0].A.el[0][i]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].a[0].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->C[n].b[0]      .el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].a[0].B    .el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].b[0].el[0].el[0][i]);
        }

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[0].C.el[1][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->B[n].a[0].el[0].el[1][i]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[0].D    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[0].el[1].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].a[0].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->B[n].b[0]      .el[1][i]);
        }

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->C[n].a[0].el[0].el[1][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[np1].a[0].C.el[0][i]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].a[0].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->C[n].b[0]      .el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].a[0].D    .el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].b[0].el[1].el[0][i]);
        }
    }
else  /* f == 1, yl_xg */
    {
    nm1 = n-1;
    
    /*  W equations  */
    for (i = 0; i < 2; i++)
        {
	fscanf(fp, "%lf", &M->AF[0].c[nm1].el[i]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->AF[0].a[n].el[i][j]);
	fscanf(fp, "%lf", &M->AF[0].b[n].el[i]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->BF[0].c[nm1].el[i][j]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M->BF[0].a[n].el[j].el[i][k]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->BF[0].b[n].el[i][j]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M->CF[0].c[nm1].el[0]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->CF[0].a[n].el[0][j]);
    fscanf(fp, "%lf", &M->CF[0].b[n].el[0]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].c[nm1].el[0].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].A.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].B.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].b[n].el[0].el[0][j]);

    /*  NE equation  */
    fscanf(fp, "%lf", &M->CF[0].c[nm1].el[1]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->CF[0].a[n].el[1][j]);
    fscanf(fp, "%lf", &M->CF[0].b[n].el[1]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].c[nm1].el[1].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].C.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].D.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].b[n].el[1].el[0][j]);
    }
}  /* end of fucntion load_W_el */


/**********************************************************************
 *  FUNCTION:  load_int_el -- loads an interior element of a collmat  *
 *  M = collmat                                                       *
 *  v = vector from which the collmat is loaded                       *
 *  f = 0, use xl_yg numbering                                        *
 *  f = 1, use yl_xg numbering                                        *
 *  n, m = indices for which element                                  *
 **********************************************************************/
void load_int_el (M, fp, f, n, m)
struct collmat *M;
FILE *fp;
int f, n, m;
{
int i, j, np1, mm1, mp1, nm1;

if (f == 0)  /* xl_yg */
    {
    np1 = n+1;
    mm1 = m-1;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[0].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .A    .el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->B[n].a[m]  .el[0].el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[m].B    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[m].el[0].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].a[m].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->B[n].b[m]      .el[0][i]);
        }

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->C[n].a[m]  .el[0].el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].c[mm1].el[0].el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].a[m]  .A    .el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].a[m].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->C[n].b[m]      .el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].a[m].B    .el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].b[m].el[0].el[0][i]);
        }

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .C    .el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->B[n].a[m]  .el[0].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[m].D    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[m].el[1].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].a[m].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->B[n].b[m]      .el[1][i]);
        }

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->C[n].a[m]  .el[0].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].c[mm1].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].a[m]  .C    .el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].a[m].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->C[n].b[m]      .el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].a[m].D    .el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].b[m].el[1].el[0][i]);
        }
    }
else  /* f == 1, yl_xg */
    {
    mp1 = m+1;
    nm1 = n-1;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[0].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].A.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].B.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[n].el[0].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].c[nm1].el[0][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->B[m].a[n].el[i].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].b[n].el[0][j]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[1].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].C.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].D.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[n].el[1].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].c[nm1].el[1][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->B[m].a[n].el[i].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].b[n].el[1][j]);

    /*  SE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].c[nm1].el[0][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->C[m].a[n].el[i].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].b[n].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].c[nm1].el[0].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].A.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].B.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].b[n].el[0].el[0][j]);

    /*  NE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].c[nm1].el[1][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->C[m].a[n].el[i].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].b[n].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].c[nm1].el[1].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].C.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].D.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].b[n].el[1].el[0][j]);
    }
}  /* end of function load_int_el */


/***************************************************************
 *  FUNCTION:  load_E_el -- loads the E elements of a collmat  *
 *  M = collmat                                                *
 *  v = vector from which the collmat is loaded                *
 *  f = 0, use xl_yg numbering                                 *
 *  f = 1, use yl_xg numbering                                 *
 *  n, m = indices for which E element                         *
 ***************************************************************/
void load_E_el (M, fp, f, n, m)
struct collmat *M;
FILE *fp;
int f, n, m;
{
int i, j, k, np1, mm1, nm1;

if (f == 0)  /* xl_yg */
    {
    np1 = n+1;
    mm1 = m-1;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[0].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .A    .el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->B[n].a[m]  .el[0].el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[m].B    .el[1][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->B[n].a[m].el[1].el[0][i]);
    
    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->C[n].a[m]  .el[0].el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].c[mm1].el[0].el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].a[m]  .A    .el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->C[n].a[m].el[1].el[0][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[np1].a[m].B    .el[0][i]);

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .C    .el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->B[n].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->B[n].a[m]  .el[0].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[m].D    .el[1][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->B[n].a[m].el[1].el[1][i]);

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->C[n].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->C[n].a[m]  .el[0].el[1][i]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[np1].c[mm1].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->A[np1].a[m]  .C    .el[0][i]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->C[n].a[m].el[1].el[1][i]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[np1].a[m].D    .el[0][i]);
    }
else  /* f == 1, yl_xg */
    {
    nm1 = n-1;
    
    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[0].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].A.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].B.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[n].el[0].el[1][j]);
    fscanf(fp, "%lf", &M->BL[0].c[nm1].el[0]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->BL[0].a[n].el[0][j]);
    fscanf(fp, "%lf", &M->BL[0].b[n].el[0]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[1].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].C.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].D.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].b[n].el[1].el[1][j]);
    fscanf(fp, "%lf", &M->BL[0].c[nm1].el[1]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->BL[0].a[n].el[1][j]);
    fscanf(fp, "%lf", &M->BL[0].b[n].el[1]);
    
    /*  E equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->CL[0].c[nm1].el[i][j]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M->CL[0].a[n].el[j].el[i][k]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->CL[0].b[n].el[i][j]);
	fscanf(fp, "%lf", &M->AL[0].c[nm1].el[i]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->AL[0].a[n].el[i][j]);
	fscanf(fp, "%lf", &M->AL[0].b[n].el[i]);
        }
    }
}  /* end of function load_E_el */


/****************************************************************
 *  FUNCTION:  load_NW_el -- loads the NW element of a collmat  *
 *  M = collmat                                                 *
 *  v = vector from which the collmat is loaded                 *
 *  f = 0, use xl_yg numbering                                  *
 *  f = 1, use yl_xg numbering                                  *
 *  n = index for which E element                               *
 ***************************************************************/
void load_NW_el (M, fp, f, n)
struct collmat *M;
FILE *fp;
int f, n;
{
int i, j, k, nm1;

if (f == 0)  /* xl_yg */
    {
    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[0].A.el[1][i]);
    fscanf(fp, "%lf", &M->BL[0].a[0].el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[0].B    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[0].el[0].el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].a[0].el[0][1]);
    fscanf(fp, "%lf", &M->BL[0].b[0].el[0]);

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->CL[0].a[0].el[0].el[0][i]);
    fscanf(fp, "%lf", &M->AL[0].a[0].el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].a[0].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->CL[0].b[0]      .el[0][i]);
	}
    fscanf(fp, "%lf", &M->AL[0].a[0].el[0][1]);
    fscanf(fp, "%lf", &M->AL[0].b[0].el[0]);

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[0].C.el[1][i]);
    fscanf(fp, "%lf", &M->BL[0].a[0].el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[0].D    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[0].el[1].el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].a[0].el[1][1]);
    fscanf(fp, "%lf", &M->BL[0].b[0].el[1]);
    
    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->CL[0].a[0].el[0].el[1][i]);
    fscanf(fp, "%lf", &M->AL[0].a[0].el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].a[0].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->CL[0].b[0]      .el[1][i]);
        }
    fscanf(fp, "%lf", &M->AL[0].a[0].el[1][1]);
    fscanf(fp, "%lf", &M->AL[0].b[0].el[1]);
    }
else  /* f == 1, yl_xg */
    {
    nm1 = n-1;

    /*  W equations  */
    for (i = 0; i < 2; i++)
        {
	fscanf(fp, "%lf", &M->AF[0].c[nm1].el[i]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->AF[0].a[n].el[i][j]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->BF[0].c[nm1].el[i][j]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M->BF[0].a[n].el[j].el[i][k]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M->CF[0].c[nm1].el[0]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->CF[0].a[n].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].c[nm1].el[0].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].A.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].B.el[0][j]);

    /*  NE equation  */
    fscanf(fp, "%lf", &M->CF[0].c[nm1].el[1]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->CF[0].a[n].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].c[nm1].el[1].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].C.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[0].a[n].D.el[0][j]);
    }
}  /* end of function load_NW_el */


/***************************************************************
 *  FUNCTION:  load_N_el -- loads the N elements of a collmat  *
 *  M = collmat                                                *
 *  v = vector from which the collmat is loaded                *
 *  f = 0, use xl_yg numbering                                 *
 *  f = 1, use yl_xg numbering                                 *
 *  n, m = indices for which N element                         *
 ***************************************************************/
void load_N_el (M, fp, f, n, m)
struct collmat *M;
FILE *fp;
int f, n, m;
{
int i, j, mm1, mp1, nm1;

if (f == 0)  /* xl_yg */
    {
    mm1 = m-1;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[0].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .A    .el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->BL[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[m].B    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[m].el[0].el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].a[m].el[0][1]);
    fscanf(fp, "%lf", &M->BL[0].b[m].el[0]);

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->CL[0].a[m]  .el[0].el[0][i]);
        }
    fscanf(fp, "%lf", &M->AL[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->AL[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].a[m].el[1].el[0][i]);
	fscanf(fp, "%lf", &M->CL[0].b[m]      .el[0][i]);
        }
    fscanf(fp, "%lf", &M->AL[0].a[m].el[0][1]);
    fscanf(fp, "%lf", &M->AL[0].b[m].el[0]);

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .C    .el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->BL[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].a[m].D    .el[1][i]);
	fscanf(fp, "%lf", &M->A[n].b[m].el[1].el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].a[m].el[1][1]);
    fscanf(fp, "%lf", &M->BL[0].b[m].el[1]);

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->CL[0].a[m]  .el[0].el[1][i]);
        }
    fscanf(fp, "%lf", &M->AL[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->AL[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].a[m].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->CL[0].b[m]      .el[1][i]);
        }
    fscanf(fp, "%lf", &M->AL[0].a[m].el[1][1]);
    fscanf(fp, "%lf", &M->AL[0].b[m].el[1]);
    }
else  /* f == 1, yl_xg */
    {
    mp1 = m+1;
    nm1 = n-1;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[0].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].A.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].B.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].c[nm1].el[0][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->B[m].a[n].el[i].el[0][j]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[1].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].C.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].D.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->B[m].c[nm1].el[1][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->B[m].a[n].el[i].el[1][j]);

    /*  SE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].c[nm1].el[0][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->C[m].a[n].el[i].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].c[nm1].el[0].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].A.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].B.el[0][j]);

    /*  NE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->C[m].c[nm1].el[1][j]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->C[m].a[n].el[i].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].c[nm1].el[1].el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].C.el[0][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[mp1].a[n].D.el[0][j]);
    }
}  /*  end of function load_N_el  */


/****************************************************************
 *  FUNCTION:  load_NE_el -- loads the NE element of a collmat  *
 *  M = collmat                                                 *
 *  v = vector from which the collmat is loaded                 *
 *  f = 0, use xl_yg numbering                                  *
 *  f = 1, use yl_xg numbering                                  *
 *  n, m = indices for which element                            *
 ****************************************************************/
void load_NE_el (M, fp, f, n, m)
struct collmat *M;
FILE *fp;
int f, n, m;
{
int i, j, k, mm1, nm1;

if (f == 0)  /* xl_yg */
    {
    mm1 = m-1;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[0].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .A    .el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->BL[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[m].B    .el[1][i]);
    fscanf(fp, "%lf", &M->BL[0].a[m].el[0][1]);

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].c[mm1]      .el[0][i]);
	fscanf(fp, "%lf", &M->CL[0].a[m]  .el[0].el[0][i]);
        }
    fscanf(fp, "%lf", &M->AL[0].c[mm1].el[0]);
    fscanf(fp, "%lf", &M->AL[0].a[m]  .el[0][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->CL[0].a[m].el[1].el[0][i]);
    fscanf(fp, "%lf", &M->AL[0].a[m].el[0][1]);

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->A[n].c[mm1].el[1].el[1][i]);
	fscanf(fp, "%lf", &M->A[n].a[m]  .C    .el[1][i]);
        }
    fscanf(fp, "%lf", &M->BL[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->BL[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->A[n].a[m].D    .el[1][i]);
    fscanf(fp, "%lf", &M->BL[0].a[m].el[1][1]);

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M->CL[0].c[mm1]      .el[1][i]);
	fscanf(fp, "%lf", &M->CL[0].a[m]  .el[0].el[1][i]);
        }
    fscanf(fp, "%lf", &M->AL[0].c[mm1].el[1]);
    fscanf(fp, "%lf", &M->AL[0].a[m]  .el[1][0]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M->CL[0].a[m].el[1].el[1][i]);
    fscanf(fp, "%lf", &M->AL[0].a[m].el[1][1]);
    }
else  /* f == 1, yl_xg */
    {
    nm1 = n-1;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[0].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].A.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].B.el[1][j]);
    fscanf(fp, "%lf", &M->BL[0].c[nm1].el[0]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->BL[0].a[n].el[0][j]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].c[nm1].el[1].el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].C.el[1][j]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->A[m].a[n].D.el[1][j]);
    fscanf(fp, "%lf", &M->BL[0].c[nm1].el[1]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M->BL[0].a[n].el[1][j]);

    /*  E equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->CL[0].c[nm1].el[i][j]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M->CL[0].a[n].el[j].el[i][k]);
	fscanf(fp, "%lf", &M->AL[0].c[nm1].el[i]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M->AL[0].a[n].el[i][j]);
        }
    }
}  /* end of function load_NE_el */


/*****************************
 *  FUNCTION:  file2collmat  *
 *****************************/
void file2collmat (M, filename, xm, ym, flag, loc, glob)
struct collmat *M;
char filename[80]; 
int xm, ym, flag, loc, glob;
{
int m, n,
    globm2 = glob-2,
    globm3 = glob-3,
    locm2  = loc-2,
    locm1  = loc-1;

FILE *fp;
  
fp = fopen(filename, "r");

/*
if (xm < ym)
    {
    if (flag == 0)   xl_yg 
        {
	load_SW_el (M, fp, flag);
	for (m = 1; m <= locm2; m++)
	    load_S_el (M, fp, flag, m);
	load_SE_el (M, fp, flag, locm1);
	for (n = 0; n <= globm3; n++)
	    {
	    load_W_el (M, fp, flag, n);
	    for (m = 1; m <= locm2; m++)
	        load_int_el (M, fp, flag, n, m);
	    load_E_el (M, fp, flag, n, locm1);
	    }
	load_NW_el (M, fp, flag, globm2);
	for (m = 1; m <= locm2; m++)
	    load_N_el (M, fp, flag, globm2, m);
	load_NE_el (M, fp, flag, globm2, locm1);
        }

    else   flag == 1, yl_xg 
        {
	load_SW_el (M, fp, flag);
	for (m = 0; m <= globm3; m++)
	    load_S_el (M, fp, flag, m);
	load_SE_el (M, fp, flag, globm2);
	for (n = 1; n <= locm2; n++)
	    {
	    load_W_el (M, fp, flag, n);
	    for (m = 0; m <= globm3; m++)
	        load_int_el (M, fp, flag, n, m);
 	    load_E_el (M, fp, flag, n, globm2);
	    }
	load_NW_el (M, fp, flag, locm1);
	for (m = 0; m <= globm3; m++)
	    load_N_el (M, fp, flag, locm1, m);
	load_NE_el (M, fp, flag, locm1, globm2);
        }
    }

else   xm >= ym 
    {
    if (flag == 0)   xl_yg 
        {
*/
	load_SW_el (M, fp, flag);
	for (n = 0; n <= globm3; n++)
	    load_W_el (M, fp, flag, n);
	load_NW_el (M, fp, flag, globm2);
	for (m = 1; m <= locm2; m++)
	    {
	    load_S_el (M, fp, flag, m);
	    for (n = 0; n <= globm3; n++)
	        load_int_el (M, fp, flag, n, m);
	    load_N_el (M, fp, flag, globm2, m);
	    }
	load_SE_el (M, fp, flag, locm1);
	for (n = 0; n <= globm3; n++)
	    load_E_el (M, fp, flag, n, locm1);
	load_NE_el (M, fp, flag, globm2, locm1);
/*
        }

    else   flag == 1, yl_xg 
        {
	load_SW_el (M, fp, flag);
	for (n = 1; n <= locm2; n++)
	    load_W_el (M, fp, flag, n);
	load_NW_el (M, fp, flag, locm1);
	for (m = 0; m <= globm3; m++)
	    {
	    load_S_el (M, fp, flag, m);
	    for (n = 1; n <= locm2; n++)
	        load_int_el (M, fp, flag, n, m);
	    load_N_el (M, fp, flag, locm1, m);
	    }
	load_SE_el (M, fp, flag, globm2);
	for (n = 1; n <= locm2; n++)
	    load_E_el (M, fp, flag, n, globm2);
	load_NE_el (M, fp, flag, locm1, globm2);	
        }
    }
*/
fclose(fp);
}  /* end of function file2collmat */


/********************************************************
 *  FUNCTION:  mult_collvec_collmat -- computes R = MU  *
 *  INPUTS: M = known collmat                           *
 *          U = known vector                            *
 *          R = result of M*U                           *
 *          flag: flag = 0 use M; flag = 1 use M'       *
 ********************************************************/
void mult_collvec_collmat (R, M, U, xm, ym)
struct collvec *R, *U;
struct collmat *M;
int xm, ym;
{
struct kvec2 tempk1, tempk2;
struct kvec4 k41, k42, k43;
int i, im1;

dkvec2(&tempk1, xm);
dkvec2(&tempk2, xm);
dkvec4(&k41, xm);
dkvec4(&k42, xm);
dkvec4(&k43, xm);

/********************************
 *  compute RF = AF*UF + BF*U0  *
 ********************************/
mult_M24(&tempk1, &M->BF[0], &U->M[0], xm);
mult_tri2xu(&tempk2, &M->AF[0], &U->F[0], xm);
addkv2(&R->F[0], &tempk1, &tempk2, xm);

/****************************************
 *  compute R0 = CF*UF + A0*U0 + B0*U1  *
 ****************************************/
mult_M42(&tempk1, &M->CF[0], &U->F[0], xm);
mult_M24(&tempk2, &M->B[0], &U->M[1], xm);
seam(&k41, &tempk1, &tempk2, xm);
mult_M44(&k42, &M->A[0], &U->M[0], xm);
addkv4(&R->M[0], &k41, &k42, xm);

/************************************************************
 *  compute R(i) = C(i-1)*U(i-1) + A(i)*U(i) + B(i)*U(i+1)  *
 ************************************************************/
for (i = 1; i <= ym-3; i++)
    {
    im1 = i-1;
    mult_M24(&tempk1, &M->C[im1], &U->M[im1], xm);
    mult_M24(&tempk2, &M->B[i], &U->M[i+1], xm);
    seam(&k41, &tempk1, &tempk2, xm);
    mult_M44(&k42, &M->A[i], &U->M[i], xm);
    addkv4(&R->M[i], &k41, &k42, xm);
    }

/*********************************************************************
 *  compute R(yn-3) = C(yn-4)*U(yn-4) + A(yn-3)*U(yn-3) + B(L)*U(L)  *
 *********************************************************************/
i = ym-2;
im1 = i-1;
mult_M24(&tempk1, &M->C[im1], &U->M[im1], xm);
mult_M42(&tempk2, &M->BL[0], &U->L[0], xm);
seam(&k41, &tempk1, &tempk2, xm);
mult_M44(&k42, &M->A[i], &U->M[i], xm);
addkv4(&R->M[i], &k41, &k42, xm);

/*********************************************
 *  compute R(L) = C(L)*U(yn-3) + A(L)*U(L)  *
 *********************************************/
mult_M24(&tempk1, &M->CL[0], &U->M[ym-2], xm);
mult_tri2xu(&tempk2, &M->AL[0], &U->L[0], xm);
addkv2(&R->L[0], &tempk1, &tempk2, xm);

free_kvec2(&tempk1);
free_kvec2(&tempk2);
free_kvec4(&k41);
free_kvec4(&k42);
free_kvec4(&k43);

}  /* end of function mult_collvec_collmat */


/****************************************************************
 *  FUNCTION:  copy_collvec -- copies one collvec into another  *
 ****************************************************************/
void copy_collvec (A, B, xm, ym)
struct collvec *A, *B;
int xm, ym;
{
int i, j, k, l;

for (j = 0; j < xm; j++)
    for (k = 0; k < 2; k++)
        {
	A->F[0].el[j].el[k] = B->F[0].el[j].el[k];
	A->L[0].el[j].el[k] = B->L[0].el[j].el[k];
        }

for (i = 0; i <= ym-2; i++)
    for (j = 0; j < xm; j++)
        for (k = 0; k < 2; k++)
	    for (l = 0; l < 2; l++)
	        A->M[i].el[j].el[k].el[l] = B->M[i].el[j].el[k].el[l];
}  /* end of function copy_collvec */


/****************************************************
 *  FUNCTION: sub_collvec -- subtract two collvecs  *
 *  INPUTS:  r, a, b   (r = a-b)                    *
 ****************************************************/
void sub_collvec (R, A, B, xm, ym)
struct collvec *R, *A, *B;
int xm, ym;
{
int i;

subkv2(&R->F[0], &A->F[0], &B->F[0], xm);
for (i = 0; i < ym-1; i++)
    subkv4(&R->M[i], &A->M[i], &B->M[i], xm);
subkv2(&R->L[0], &A->L[0], &B->L[0], xm);
}  /* end of function sub_collvec */


/****************************************************
 *  FUNCTION: psolve -- solves matrix equation      *
 *         A = D-L (i.e., A = Gauss-Seidel matrix)  *
 ****************************************************/
void psolve (x, A, b, xm, ym)
struct collvec *x, *b;
struct collmat *A;
int xm, ym;
{
struct kvec2 tempk1, tempk2;
struct kvec4 k41, k42, k43;
int i, j, k, l, m, im1, xmm1=xm-1;

dkvec2(&tempk1, xm);
dkvec2(&tempk2, xm);
dkvec4(&k41, xm);
dkvec4(&k42, xm);
dkvec4(&k43, xm);

if ((ym-1) % 2 == 0)   /* yn even */
    {
    /*********
     *  RED  *
     *********/
    tridiag2(&x->F[0], &A->AF[0], &b->F[0], xm);

    for (i = 1; i <= ym-2; i += 2)
        tridiag4(&x->M[i], &A->A[i], &b->M[i], xm);
	
    /***********
     *  BLACK  *
     ***********/
    for (i = 0; i <= ym-3; i += 2)
        {
	if (i == 0)
	    {
	    mult_M42(&tempk1, &A->CF[0], &x->F[0], xm);
	    mult_M24(&tempk2, &A->B[0], &x->M[1], xm);
	    seam(&k41, &tempk1, &tempk2, xm);
	    subkv4(&k42, &b->M[0], &k41, xm);
	    }
	else
	    {
	    im1 = i-1;
	    mult_M24(&tempk1, &A->C[im1], &x->M[im1], xm);
	    mult_M24(&tempk2, &A->B[i], &x->M[i+1], xm);
	    seam(&k41, &tempk1, &tempk2, xm);
	    subkv4(&k42, &b->M[i], &k41, xm);
	    }
	    
	tridiag4(&x->M[i], &A->A[i], &k42, xm);
        }

    mult_M24(&tempk1, &A->CL[0], &x->M[ym-2], xm);
    subkv2(&tempk2, &b->L[0], &tempk1, xm);
    tridiag2(&x->L[0], &A->AL[0], &tempk2, xm);
    }
else /* yn odd */
    {
    /*********
     *  RED  *
     *********/
    tridiag2(&x->F[0], &A->AF[0], &b->F[0], xm);
	
    for (i = 1; i <= ym-3; i += 2)
        tridiag4(&x->M[i], &A->A[i], &b->M[i], xm);
	
    tridiag2(&x->L[0], &A->AL[0], &b->L[0], xm);
	
    /***********
     *  BLACK  *
     ***********/
    for (i = 0; i <= ym-2; i += 2)
        {
	if (i == 0)
	    {
	    mult_M42(&tempk1, &A->CF[0], &x->F[0], xm);
	    mult_M24(&tempk2, &A->B[0], &x->M[1], xm);
	    seam(&k41, &tempk1, &tempk2, xm);
	    subkv4(&k42, &b->M[0], &k41, xm);
	    }
	else if (i == ym-2)
	    {
	    im1 = i-1;
	    mult_M24(&tempk1, &A->C[im1], &x->M[im1], xm);
	    mult_M42(&tempk2, &A->BL[0], &x->L[0], xm);
	    seam(&k41, &tempk1, &tempk2, xm);
	    subkv4(&k42, &b->M[i], &k41, xm);
	    }
	else
	    {
	    im1 = i-1;
	    mult_M24(&tempk1, &A->C[im1], &x->M[im1], xm);
	    mult_M24(&tempk2, &A->B[i], &x->M[i+1], xm);
	    seam(&k41, &tempk1, &tempk2, xm);
	    subkv4(&k42, &b->M[i], &k41, xm);
	    }
	    
	tridiag4(&x->M[i], &A->A[i], &k42, xm);
        }
    }

free_kvec2(&tempk1);
free_kvec2(&tempk2);
free_kvec4(&k41);
free_kvec4(&k42);
free_kvec4(&k43);
}  /* end of function psolve */


/***********************************************
 *  FUNCTION:  mult_collvec_collmat_trick --   *
 *                 computes T = MZ  (using S)  *
 ***********************************************/
void mult_collvec_collmat_trick (T, M, Z, S, xm, ym)
struct collvec *T, *Z, *S;
struct collmat *M;
int xm, ym;
{
int i, im1, ymm2 = ym-2;
struct kvec2 tempk1, tempk2;

dkvec2(&tempk1, xm);
dkvec2(&tempk2, xm);

if (ym % 2 == 0)  /* ym even */
    {
    /***********
     *  BLACK  *
     ***********/
    for (i = 0; i <= ymm2; i+=2)
        copy_kvec4(&T->M[i], &S->M[i], xm);

    /*********
     *  RED  *
     *********/
    mult_M24(&T->F[0], &M->BF[0], &Z->M[0], xm);
    addkv2(&T->F[0], &T->F[0], &S->F[0], xm);

    for (i = 1; i <= ym-3; i+=2)
        {
	im1 = i - 1;
	mult_M24(&tempk1, &M->C[im1], &Z->M[im1], xm);
	mult_M24(&tempk2, &M->B[i],   &Z->M[i+1], xm);
	seam(&T->M[i], &tempk1, &tempk2, xm);
	addkv4(&T->M[i], &T->M[i], &S->M[i], xm);
        }

    mult_M24(&T->L[0], &M->CL[0], &Z->M[ymm2], xm);
    addkv2(&T->L[0], &T->L[0], &S->L[0], xm);
    }
else  /* ym odd */
    {
    printf("ym must be even\nExiting from function mult_collvec_collmat_trick\n");
    exit(0);
    }

free_kvec2(&tempk1);
free_kvec2(&tempk2);

}  /* end of function mult_collvec_collmat_trick */


/***************************************************************
 *  FUNCTION: pbcgstab -- solves Ax=b by preconditioned        *
 *                   bi-conjugate gradient method stabilized   *
 ***************************************************************/
void pbcgstab (x, A, b, xm, ym, maxiter, tol, iter)
struct collvec *x, *b;
struct collmat *A;
int xm, ym, maxiter, iter;
double tol;
{
struct timeval t_s, t_f;
struct collvec r, rbar, v, p, y, s, t, z;
double rho, oldrho, alpha, omega, beta, reltol, nm, secs;

dcollvec(&r,    xm, ym);
dcollvec(&rbar, xm, ym);
dcollvec(&v,    xm, ym);
dcollvec(&p,    xm, ym);
dcollvec(&y,    xm, ym);
dcollvec(&s,    xm, ym);
dcollvec(&t,    xm, ym);
dcollvec(&z,    xm, ym);

mult_collvec_collmat(&r, A, x, xm, ym);
sub_collvec(&r, b, &r, xm, ym);
copy_collvec(&rbar, &r, xm, ym);
reltol = tol * tol * dot_collvec(&r, &r, xm, ym);

printf("\nreltol = %e\n\n", reltol);

rho = 1;
alpha = 1;
omega = 1;

iter = 0;

gettimeofday(&t_s, (struct timezone*)0);

while (iter < maxiter)
    {
    iter++;

    oldrho = rho;
    rho = dot_collvec(&rbar, &r, xm, ym);
    beta = (rho/oldrho) * (alpha/omega);
    
    saxpy_collvec(&p, -omega, &v, &p, xm, ym);
    saxpy_collvec(&p, beta,   &p, &r, xm, ym);
    
    psolve(&y, A, &p, xm, ym);    
    mult_collvec_collmat_trick(&v, A, &y, &p, xm, ym);

    alpha = rho / dot_collvec(&rbar, &v, xm, ym);
    saxpy_collvec(&s, -alpha, &v, &r, xm, ym);

    psolve(&z, A, &s, xm, ym);    
    mult_collvec_collmat_trick(&t, A, &z, &s, xm, ym);

    omega = dot_collvec(&t, &s, xm, ym) / dot_collvec(&t, &t, xm, ym);
    
    saxpy_collvec(x,   alpha, &y, x,  xm, ym);
    saxpy_collvec(x,   omega, &z, x,  xm, ym);
    saxpy_collvec(&r, -omega, &t, &s, xm, ym);
    
    nm = dot_collvec(&r, &r, xm, ym);

    printf("%ld   %e\n", iter, nm);

    if (nm <= reltol)
        break;
    }

gettimeofday(&t_f, (struct timezone*)0);
secs = ((t_f.tv_sec-t_s.tv_sec)*1000000 + t_f.tv_usec - t_s.tv_usec)/1000000.;
printf("time for main loop = %lf seconds for single processor\n", secs);

if (iter >= maxiter)
    {
    printf("\n\007Maximum number of %d iterations exceeded.\n", maxiter);
    printf("Exiting program...\n\n");
    exit(0);
    }

free_collvec(&r   , ym);
free_collvec(&rbar, ym);
free_collvec(&v   , ym);
free_collvec(&p   , ym);
free_collvec(&y   , ym);
free_collvec(&s   , ym);
free_collvec(&t   , ym);
free_collvec(&z   , ym);
}  /* end of function pbcgstab */


/******************************************************
 *  cv2vA -- one way to load a collvec into a vector  *
 *  INPUTS:  v = vector                               *
 *           K = collvec                              *
 *           loc, glob = number of elements in        *
 *                       coordinate directions        *
 ******************************************************/
void cv2vA (v, K, loc, glob)
double *v;
struct collvec *K;
int loc, glob;
{
int i, j, k, l, c=0;

for (i = 0; i < loc; i++)
    for (j = 0; j < 2; j++)
        v[c++] = K->F[0].el[i].el[j];

for (i = 0; i <= glob-2; i++)
    for (l = 0; l < 2; l++)
        for (j = 0; j < loc; j++)
	    for (k = 0; k < 2; k++)
	        v[c++] = K->M[i].el[j].el[k].el[l];

for (i = 0; i < loc; i++)
    for (j = 0; j < 2; j++)
        v[c++] = K->L[0].el[i].el[j];
}  /* end of function cv2vA */


/************************************************************
 *  cv2vB -- the other way to load a collvec into a vector  *
 *  INPUTS:  v = vector                                     *
 *           K = collvec                                    *
 *           loc, glob = number of elements in              *
 *                       coordinate directions              *
 ************************************************************/
void cv2vB (v, K, loc, glob)
double *v;
struct collvec *K;
int loc, glob;
{
int i, j, k, l, c=0;

for (l = 0; l < loc; l++)
    for (i = 0; i < 2; i++)
        {
	v[c++] = K->F[0].el[l].el[i];
	for (j = 0; j <= glob-2; j++)
	    for (k = 0; k < 2; k++)
	        v[c++] = K->M[j].el[l].el[i].el[k];
	v[c++] = K->L[0].el[l].el[i];
        }
}  /* end of fucntion cv2vB */


/********************************************************
 *  FUNCTION:  collvec2vec -- given a collvec           *
 *             load this collvec into a vector          *
 *             The numbering is such that               *
 *                                                      *
 *  INPUTS: v = vector to be loaded                     *
 *          K = given collvec                           *
 *                                                      *
 *          xm = number of elements in x-direction      *
 *          ym = number of elements in y-direction      *
 ********************************************************/
void collvec2vec (v, K, loc, glob, flag, xmm, ymm)
struct collvec *K;
double *v;
int loc, glob, xmm, ymm;
long flag;
{
if (xmm < ymm)
    {
    if (flag == 0)  /* xl_yg */
        cv2vA (v, K, loc, glob);
    else  /* flag == 1, yl_xg */
        cv2vB (v, K, loc, glob);
    }
else  /* xmm >= ymm */
    {
    if (flag == 1)  /* yl_xg */
        cv2vA (v, K, loc, glob);
    else  /* flag == 0, xl_yg */
        cv2vB (v, K, loc, glob);
    }
}  /* end of fucntion collvec2vec */


/**************************
 *  FUNCTION:  run_1proc  *
 **************************/
void run_1proc(solvec, xm, ym, num, nproc, maxiter, eps, Afile, xfile, bfile)
int xm, ym, num, nproc, maxiter;
double *solvec, eps;
char *Afile, *xfile, *bfile;
{
int loc, glob, iter;

struct collvec U, K;
struct collmat L;

/*********************
 * define loc, glob  *
 *********************/
if (num == 0)
    {
    loc  = xm;
    glob = ym;
    }
else
    {
    loc  = ym;
    glob = xm;
    }
  
/******************
 *  read in data  *
 ******************/
dcollvec(&U, loc, glob);
dcollvec(&K, loc, glob);
dcollmat(&L, loc, glob);
file2collvec(&U, xfile, xm, ym, num, loc, glob);
file2collvec(&K, bfile, xm, ym, num, loc, glob);
file2collmat(&L, Afile, xm, ym, num, loc, glob);

/***********
 *  solve  *
 ***********/
iter = 0;
pbcgstab(&U, &L, &K, loc, glob, maxiter, eps, iter);

collvec2vec (solvec, &K, loc, glob, num, xm, ym);

/*****************
 *  de-allocate  *
 *****************/
free_collvec(&U, glob);
free_collvec(&K, glob);
free_collmat(&L, glob);

}  /* end of function run_1proc */