#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "include.h"
#include "spectrum.h"  /* the prototype for freeqstore in here */

/* Memory used by iovect and laso */

static struct {
  doublereal *real;  /* real and complex point to the same memory */
  doublecomplex *complex;
  size_t length; /* length of array in bytes */
} qstore={NULL,NULL,0};

/*
   This is the routine used by lanso* to store Lanczos vectors
*/

extern fsub STORE(integer *n, integer *isw, integer *j, doublereal *s){

  size_t trylength=(*j)*(*n)*sizeof(*s);
  if(trylength>qstore.length){
    qstore.real=realloc(qstore.real,qstore.length=trylength);
    qstore.complex=(doublecomplex *) qstore.real;
     if(qstore.real == NULL){
       fprintf(stderr,__FILE__":  Memory error in STORE\n");
      exit(55);
    }
  }

#if 0 /* debug: check orthogonality against previous vectors */
  if(*isw==1){
    int i,k,l=0;
    doublereal x,*p;
    for(i=1; i<*j; i++){
      x=0.0;
      p=qstore.real+(i-1)*(*n);
      for(k=0; k<*n; k++)x+=p[k]*s[k];
      /* Wu & Simon, p. 9, this is "partial orthogonality"  */
      if(fabs(x) > sqrt(DBL_EPSILON)){  
	printf("   Loss of orthogonality, vectors J=%li and %i:",*j,i);
	printf("  %g\n",x);
	if(l++>3){ printf("   ...\n"); break; }
      }
    }
  }
#endif

  if (*isw == 1) {
    memcpy(qstore.real+((*j)-1)*(*n),s,*n*sizeof(*s));
  } else if (*isw == 2) {
    memcpy(s,qstore.real+((*j)-1)*(*n),*n*sizeof(*s));
  } else {
    fprintf(stderr,__FILE__":  error in zSTORE, isw=%li\n",*isw);
    exit(55);
  }
  return SUBRETURN;
}

/* complex version of the above */

extern fsub ZSTORE(integer *n, integer *isw, integer *j, doublecomplex *s){

  size_t trylength=sizeof(*s)*(*j)*(*n);
  if(trylength>qstore.length){
    qstore.real=realloc(qstore.real,qstore.length=trylength);
    qstore.complex=(doublecomplex *) qstore.real;
    if(qstore.complex == NULL){
       fprintf(stderr,__FILE__":  Memory error in zSTORE, trylength=%i\n",
	       trylength);
      exit(55);
    }
  }

#if 0 /* debug: check orthogonality against previous vectors */
  if(*isw==1){
    int i,k,l=0;
    doublecomplex x,*p;
    for(i=1; i<*j; i++){
      x.re=0.0; x.im=0.0;
      p=qstore.complex+(i-1)*(*n);
      for(k=0; k<*n; k++){
	x.re+=p[k].re*s[k].re+p[k].im*s[k].im;
	x.im+=-p[k].im*s[k].re+p[k].re*s[k].im;
      }
      /* Wu & Simon, p. 9, this is "partial orthogonality"  */
      if(pow(x.re,2)+pow(x.im,2) > DBL_EPSILON){
	printf("   Loss of orthogonality, vectors J=%li and %i:",*j,i);
	printf("  %g%+g*i\n",x.re,x.im);
	if(l++>3){ printf("   ...\n"); break; }
      }
    }
  }
#endif
 
  if (*isw == 1) {
    memcpy(qstore.complex+((*j)-1)*(*n),s,*n*sizeof(*s));
  } else if (*isw == 2) {
    memcpy(s,qstore.complex+((*j)-1)*(*n),*n*sizeof(*s));
  } else {
    fprintf(stderr,__FILE__":  error in zSTORE, isw=%li\n",*isw);
    exit(55);
  }
  return SUBRETURN;
}


/*************************************************************************/

/*  This is the subroutine that stores vectors for dnlaso.
    See the documentation file "laso.doc"      */

#define DEBUGIO 0 /* flag to print out debug statements */

extern fsub iovect(integer *n, integer *m, doublereal *q, integer *j, 
		  integer *k){
  integer i, l, l1;
  size_t temp_length=sizeof(*q)* *j* *n;

  if(temp_length > qstore.length){
    qstore.real=realloc(qstore.real,qstore.length=temp_length);
    qstore.complex=(doublecomplex *) qstore.real;
    if(NULL==qstore.real){
       fprintf(stderr,__FILE__":  Memory error in STORE\n");
      exit(55);
    }
  }
  
  if(*k==0){
    for(l=0; l<*m; l++){
      l1 = *j - *m + l;
#if defined(__uxp__)
#pragma loop noalias qstore
#pragma loop disjoint qstore
#endif
      for(i=0; i<*n; i++){
	qstore.real[i+l1 * *n]=q[i+l * *n];
      }
    }
    return SUBRETURN;
  }else {
    for(l=0; l<*m; l++){
      l1 = *j - *m + l;
#if defined(__uxp__)
#pragma loop noalias q
#pragma loop disjoint q
#endif
      for (i=0; i<*n; i++){
	q[i + l * *n] = qstore.real[i + l1 * *n];
      }
    }
    return SUBRETURN;
  }
}

/* free up memory taken by qstore */

void free_qstore(void){
  free(qstore.real);
  free(qstore.complex);
  qstore.length=0;
  qstore.real=NULL;
  qstore.complex=NULL;
}