symamgsavediag.c

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


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


double SYMAMGSAVEDIAG(size_t *Fparam, integer *n, integer *ia, 
		    integer *ja, FLOAT *a)
{
  integer i,j,k,ki;
  REALS   val,ti,tmp;
  REALS   *disc;
  ILUPACKPARAM *IPparam;
  
  IPparam=(ILUPACKPARAM *)MALLOC((size_t)1*sizeof(ILUPACKPARAM),"ilupackfactor:param");
  memcpy(Fparam, &IPparam, sizeof(size_t));

  IPparam->diag=(FLOAT *)MALLOC(*n*sizeof(FLOAT),"symamgsavediag:IPparam->diag");
  disc=(REALS *)MALLOC(*n*sizeof(REALS),"symamgsavediag:disc");

  for (i=0; i<*n; i++)
      disc[i]=IPparam->diag[i]=0.0;

  ki=1;
  for (i=0; i<*n; i++) {
      for (j=ia[i]; j<ia[i+1]; j++) {
	  k=ja[j-1]-1;
	  if (k==i) {
	     IPparam->diag[i]=a[j-1];
	     if (i==0) 
		 ti=a[j-1];
	     else {
		 if (a[j-1] < ti)  
		     ti=a[j-1];
		 ki=i;
	     }
	  }
	  else {
	      val=FABS(a[j-1]);
	      disc[i]+=val;
	      // remember that only the upper triangular part is stored
	      disc[k]+=val;
	  }
      }
  }


  val=ti-disc[ki];
  for (i=1; i<*n; i++) {
      tmp=IPparam->diag[i]-disc[i]; 
      if (tmp<val) {
	  if (IPparam->diag[i]-ti < disc[ki]) 
	      val=tmp;
      }
  }
  free(disc);

  return (val);
}