dglprecsetup.c

来自「a software code for computing selected e」· C语言 代码 · 共 62 行

#include <string.h>
#include <ilupack.h>
#include <ilupackmacros.h>

#define MAX(A,B)   (((A)>=(B))?(A):(B))
#define MIN(A,B)   (((A)<=(B))?(A):(B))

void  dglprecsetup(size_t    *IPdiag, 
		   FLOAT     *a,
		   integer   *n,
		   integer   *ISEARCH,
		   REALS     *SHIFT,
		   REALS     *shift0,
		   integer   *factdgl,
		   integer   *factspd,
		   REALS     *diagmin) {

  FLOAT *diag;
  FLOAT ta,del,tmp,ati,ata;
  integer i;

  diag=(FLOAT *)MALLOC(*n*sizeof(FLOAT),"dglprecsetup:diag");
  memcpy(IPdiag, &diag, sizeof(size_t));

  *diagmin=a[0];
  ta=*diagmin;
  for (i=1; i<*n; i++) {
      *diagmin=MIN(*diagmin,a[i]);
      ta=MAX(ta,a[i]);
  }
  ati=MAX(*diagmin,-*diagmin);
  ata=MAX(ta,-ta);
  ata=MIN(ata,ati);
  del=MAX(ta-*diagmin,ata);
  
  if (*ISEARCH<=1) {
      *shift0=*diagmin-del/(*n);  
      *factspd=-1;
      if (*ISEARCH==1) {
         if (*SHIFT<*diagmin) {
	    *shift0=*SHIFT;
	    *factspd=0;    
         }
      }
      for (i=0; i<*n; i++)
        diag[i]=1.0/(a[i]-*shift0);
  }
  else {
      *factspd=0;
      del=1./del;
      for   (i=0; i<*n; i++) {
         tmp=a[i]-*shift0;
         if (tmp > 0.0)
	     tmp=MAX(tmp,del);
	 else
	     tmp=MIN(tmp,-del);
         diag[i]=1.0/tmp;
      }
  }
  *factdgl=1;
}