matrix.c

cfd求解器使用与gmsh网格的求解

    }
  }

}

void multi_prod_vector_double(int n, int Sizevec, double **Vec, double *coef, double *VecRes ){
  int k;

  zero_vector(Sizevec,VecRes);
  for(k=0;k<n;k++){
    if (coef[k]){
      prod_vector_double(Sizevec,Vec[k],coef[k]);
      add_vector_vector(Sizevec,VecRes,Vec[k]);
      prod_vector_double(Sizevec,Vec[k],1.0/coef[k]);
    }
  }
}

void multi_prod_matrix_vector(int n, int Sizevec, Matrix **Mat, double **Vec, double *VecRes ){
  int k;
  double *work;

  init_vector(Sizevec,&work);

  zero_vector(Sizevec,VecRes);
  for(k=0;k<n;k++){
    prod_matrix_vector(Mat[k],Vec[k],work);
    add_vector_vector(Sizevec,VecRes,work);
  }
}

void prodsc_vectorconj_vector (int Nb, double *U , double *V, double *proscar, double *proscai ){
  /* uconjugue * v -> proscar + i prodscai  */
  int i;
  *proscar = *proscai = 0.0 ;
  for (i=0; i<Nb; i+=2){
    *proscar += U[i] * V[i]   + U[i+1] * V[i+1] ;
    *proscai += U[i] * V[i+1] - U[i+1] * V[i]   ;
  }
}

/* ------------------------------------------------------------------------ */
/*  n o r m                                                                 */
/* ------------------------------------------------------------------------ */

void norm2_vector(int Nb, double *U, double *norm){
  prodsc_vector_vector(Nb,U,U,norm);
  *norm = sqrt(*norm);
}

void norminf_vector(int Nb, double *U, double *norm){
  int i;
  *norm = 0.;
  for(i=0; i<Nb; i++) 
    if(fabs(U[i]) > *norm) *norm = fabs(U[i]);
}

/* ------------------------------------------------------------------------ */
/*  i d e n t i t y                                                         */
/* ------------------------------------------------------------------------ */

void identity_matrix (Matrix *M){
  int i;
  zero_matrix(M);
  for (i=1;i<=M->N;i++) add_matrix_double(M,i,i,1.0);
}


/* ------------------------------------------------------------------------ */
/*  w r i t e                                                               */
/* ------------------------------------------------------------------------ */

void binary_write_matrix (Matrix *M, char *name, char *ext){
  
  int   Nb;
  FILE *pfile;
  char  filename[256];

  if(!M->N){
    Msg(WARNING, "No elements in matrix");
    return;
  }

  strcpy(filename, name);
  strcat(filename, ext);
  pfile = fopen(filename, "wb") ;

  fprintf(pfile,"%d\n",M->T);

  switch (M->T) {
  case SPARSE :
    Nb = List_Nbr(M->S.a) ;

    fprintf(pfile,"%d %d\n", M->N, Nb);

    fprintf(pfile,"%d %d %d %d %d\n", M->S.ptr->nmax, M->S.ptr->size, 
	    M->S.ptr->incr, M->S.ptr->n, M->S.ptr->isorder);
    fprintf(pfile,"%d %d %d %d %d\n", M->S.ai->nmax, M->S.ai->size, 
	    M->S.ai->incr, M->S.ai->n, M->S.ai->isorder);
    fprintf(pfile,"%d %d %d %d %d\n", M->S.jptr->nmax, M->S.jptr->size, 
	    M->S.jptr->incr, M->S.jptr->n, M->S.jptr->isorder);
    fprintf(pfile,"%d %d %d %d %d\n", M->S.a->nmax, M->S.a->size, 
	    M->S.a->incr, M->S.a->n, M->S.a->isorder);

    fwrite(M->S.ptr->array, sizeof(int), Nb, pfile);
    fwrite(M->S.ai->array, sizeof(int), Nb, pfile);
    fwrite(M->S.jptr->array, sizeof(int), M->N, pfile);
    fwrite(M->S.a->array, sizeof(double), Nb, pfile);
    break;

  case DENSE :
    fprintf(pfile,"%d\n", M->N);
    fwrite(M->F.a, sizeof(double), M->N*M->N, pfile);
    break;
  }

  fclose(pfile) ;
}


void binary_write_vector (int Nb, double *V, char *name, char *ext){
  char filename[256];
  FILE *pfile;

  strcpy(filename, name);
  strcat(filename, ext);
  pfile = fopen(filename, "wb") ;

  fwrite(V, sizeof(double), Nb, pfile);

  fclose(pfile) ;
}


void formatted_write_matrix (FILE *pfile, Matrix *M, int style){  
  int *ptr,*ai,i,j,*jptr, *ia, *ja, *ir, nnz, ierr;
  int  un=1;
  double *a;

  if(!M->N){
    Msg(WARNING, "No element in matrix");
    return;
  }

  switch (M->T) {

  case DENSE :
    if(M->notranspose)
      for(i=0 ; i<M->N ; i++)
	for(j=0 ; j<M->N ; j++)
	  fprintf(pfile,"%d %d %.16g\n", j+1, i+1, M->F.a[i*(M->N)+j]);
    else
      for(i=0 ; i<M->N ; i++)
	for(j=0 ; j<M->N ; j++)
	  fprintf(pfile,"%d %d %.16g\n", i+1, j+1, M->F.a[i*(M->N)+j]);
    break;

  case SPARSE :
    
    switch(style){
      
    case ELAP : 
      fprintf(pfile,"%d\n",M->T); 
      a = (double*)M->S.a->array;
      ai = (int*)M->S.ai->array;
      ptr = (int*)M->S.ptr->array;
      jptr = (int*)M->S.jptr->array;
      fprintf(pfile,"%d\n",M->N);
      fprintf(pfile,"%d\n",List_Nbr(M->S.a));
      for(i=0;i<M->N;i++)
	fprintf(pfile," %d",jptr[i]);
      fprintf(pfile,"\n");
      for(i=0;i<List_Nbr(M->S.a);i++)
	fprintf(pfile,"%d %d %.16g \n",ai[i],ptr[i],a[i]);
      break;
    
    case KUL :
      csr_format(&M->S, M->N);
      a  = (double*) M->S.a->array;
      ia = (int*) M->S.jptr->array;
      ja = (int*) M->S.ptr->array; 
      nnz = List_Nbr(M->S.a);
      ir = (int*) Malloc(nnz * sizeof(int));
      csrcoo_(&M->N, &un, &nnz, a, ja, ia, &nnz, a, ir, ja, &ierr);      
      for(i=0 ; i<nnz ; i++)
	fprintf(pfile,"%d  %d  %.16g\n", ir[i], ja[i], a[i]);
      restore_format(&M->S);
      break;
    
    default : 
      Msg(GERROR, "Unknown print style for formatted matrix output");
    }
    break ;
    
  default :
    Msg(GERROR, "Unknown matrix format for formatted matrix output");

  }
}


void formatted_write_vector (FILE *pfile, int Nb, double *V, int style){
  int  i;

  /* for(i=0 ; i<Nb ; i++) fprintf(pfile,"%d %.16g\n", i+1, V[i]); */
  for(i=0 ; i<Nb ; i++) fprintf(pfile,"%.16g\n", V[i]);
}


/* ------------------------------------------------------------------------ */
/*  r e a d                                                                 */
/* ------------------------------------------------------------------------ */


void binary_read_matrix (Matrix *M, char *name , char *ext){
  
  int   Nb;
  FILE *pfile;
  char  filename[256];
  
  strcpy(filename, name);
  strcat(filename, ext);
  pfile = fopen(filename, "rb") ;

  if (pfile == NULL) {
    Msg(GERROR,"Error opening file '%s'", filename);    
  }
  
  fscanf(pfile,"%d",&M->T);
  M->ILU_Exists  = 0;

  switch (M->T) {
  case SPARSE :
    fscanf(pfile,"%d %d\n", &M->N, &Nb);

    M->S.ptr  = List_Create (Nb, 1, sizeof(int));
    M->S.ai   = List_Create (Nb, 1, sizeof(int));
    M->S.jptr = List_Create (M->N, 1, sizeof(int));
    M->S.a    = List_Create (Nb, 1, sizeof(double));

    fscanf(pfile,"%d %d %d %d %d\n", &M->S.ptr->nmax, &M->S.ptr->size, 
	   &M->S.ptr->incr, &M->S.ptr->n, &M->S.ptr->isorder);
    fscanf(pfile,"%d %d %d %d %d\n", &M->S.ai->nmax, &M->S.ai->size, 
	   &M->S.ai->incr, &M->S.ai->n, &M->S.ai->isorder);
    fscanf(pfile,"%d %d %d %d %d\n", &M->S.jptr->nmax, &M->S.jptr->size, 
	   &M->S.jptr->incr, &M->S.jptr->n, &M->S.jptr->isorder);
    fscanf(pfile,"%d %d %d %d %d\n", &M->S.a->nmax, &M->S.a->size, 
	   &M->S.a->incr, &M->S.a->n, &M->S.a->isorder);

    fread(M->S.ptr->array, sizeof(int), Nb, pfile);
    fread(M->S.ai->array, sizeof(int), Nb, pfile);
    fread(M->S.jptr->array, sizeof(int), M->N, pfile);
    fread(M->S.a->array, sizeof(double), Nb, pfile);
    break;
    
  case DENSE :
    fscanf(pfile,"%d\n", &M->N);
    M->F.LU_Exist = 0;
    M->F.a  = (double*) Malloc(M->N * M->N * sizeof(double));
    M->F.lu  = (double*) Malloc(M->N * M->N * sizeof(double));
    fread(M->F.a, sizeof(double), M->N * M->N, pfile);
    break ;
  }

  fclose(pfile) ;

}


void binary_read_vector (int Nb, double **V, char *name, char *ext){
  char filename[256];
  FILE *pfile;

  strcpy(filename, name);
  strcat(filename, ext);
  pfile = fopen(filename, "rb") ;

  if (pfile == NULL) {
    Msg(GERROR, "Error opening file %s", filename);
  }

  init_vector(Nb, V);
  fread(*V, sizeof(double), Nb, pfile);

  fclose(pfile) ;
}


void formatted_read_matrix (Matrix *M, char *name , char *ext, int style){  
  int i,nnz,inb,inb2;
  double nb;
  FILE *pfile;
  char filename[256];
  
  strcpy(filename, name);
  strcat(filename, ext);
  pfile = fopen(filename, "r") ;

  if (pfile == NULL) {
    Msg(GERROR,"Error opening file  %s", filename);
  }
  
  fscanf(pfile,"%d",&M->T);
  switch (M->T) {
  case SPARSE :
    List_Reset(M->S.jptr);
    fscanf(pfile,"%d",&M->N);
    fscanf(pfile,"%d",&nnz);
    for(i=0;i<M->N;i++){
      fscanf(pfile," %d",&inb);
      List_Add(M->S.jptr,&inb);
    }
    for(i=0;i<nnz;i++){
      fscanf(pfile,"%d %d %lf \n",&inb,&inb2,&nb);
      List_Add(M->S.ai,&inb);
      List_Add(M->S.ptr,&inb2);
      List_Add(M->S.a,&nb);
    }
    
    break;
    
  case DENSE :
    fscanf(pfile,"%d",&M->N);
    for(i=0;i<(M->N)*(M->N);i++){
      fscanf(pfile,"%d %lf ", &inb, &M->F.a[i]);
    }
    break;
  }
  fclose(pfile) ;

}


void formatted_read_vector (int Nb, double *V, char *name, char *ext, int style){
  int i;
  FILE *pfile;
  char filename[256];

  strcpy(filename, name);
  strcat(filename, ext);
  pfile = fopen(filename, "r") ;

  if (pfile == NULL) {
    Msg(GERROR,"Error opening file %s", filename);
  }

  for(i=0 ; i<Nb ; i++) fscanf(pfile, "%lf", &V[i]);

  fclose(pfile) ;
}

/* ------------------------------------------------------------------------ */
/*  p r i n t _ m a t r i x _ i n f o _ X X X                               */
/* ------------------------------------------------------------------------ */

int maximum (int a, int b) {
  if (a>b) 
    return(a);
  else
    return(b);
}


void print_matrix_info_CSR (int N, int *jptr, int *ai){
  int i, j, k, l, m, n;
  
  l = n = 0;
  j = jptr[N]-1 ;
  for (i=0; i<N; i++) {
    k = jptr[i+1] - jptr[i];	  
    m = ai[jptr[i+1]-2] - ai[jptr[i]-1] + 1;
    if (l<k) l = k;
    if (n<m) n = m;
  }

  Msg(SPARSKIT, "N: %d, NZ: %d, BW max/avg: %d/%d, SW max: %d\n", 
      N, j, l, (int)(j/N), n);
}

void print_matrix_info_MSR (int N, scalar *a, int *jptr){
  int i, j, k, l, m, n;

  l = n = 0;
  j = jptr[N]-2;
  for (i=0; i<N; i++) {
    k = jptr[i+1] - jptr[i] + (a[i]?1:0) ; 
    if((jptr[i+1] - jptr[i]) == 0)
      m = (a[i]?1:0);
    else
      m = maximum ( jptr[jptr[i+1]-2] - jptr[jptr[i]-1] + 1,
		    maximum (jptr[jptr[i+1]-2] - (i+1) + 1, (i+1) - jptr[jptr[i]-1] + 1) );
    if (l<k) l = k;
    if (n<m) n = m;
  }

  Msg(SPARSKIT, "N: %d, NZ: %d, BW max/avg: %d/%d, SW max: %d\n", 
      N, j, l, (int)(j/N), n);

}

void print_matrix_info_DENSE (int N){
  Msg(SPARSKIT, "N: %d\n", N);
}


/* ------------------------------------------------------------------------ */
/*  get _ column _ in _ m a t r i x                                         */
/* ------------------------------------------------------------------------ */

void get_column_in_matrix (Matrix *M, int col, double *V){

  int     k, i, j, *ai, *jptr ;
  double  *a;
  int found;

  switch (M->T) {
  case SPARSE :
    /* csr_format transpose! 
       donc la matrice arrivant dans cette routine doit
       bel et bien etre la transposee !!! */
    if(M->changed){
      csr_format (&M->S, M->N);
      restore_format (&M->S);
      M->changed = 0 ;
    }
    jptr = (int*) M->S.jptr->array;
    a    = (double*) M->S.a->array;
    ai   = (int*) M->S.ai->array;

    for(i=0; i<M->N; i++){  /* lignes */
      found=0;
      for(k=jptr[i]-1;k<jptr[i+1]-1;k++){ /*colonne */
         if(ai[k]-1==col) { 
	   V[i]=a[k]; found=1; break; 
	 }
	 else if (ai[k]-1 > col) { 
	   break; 
	 }
       }
      if (!found) V[i]=0; 
    printf(" V[%d] = %g \n",i, V[i]); 
    } 
    break;
  case DENSE :
    if(M->notranspose){
	for(j=0; j<M->N; j++) V[j] = M->F.a[(M->N)*col+j];
    }
    else{
      for(i=0; i<M->N; i++){
	for(j=0; j<M->N; j++) V[j] = M->F.a[(M->N)*j+col];
      }
    }
    break;
  }
}



void get_element_in_matrix (Matrix *M, int row, int col, double *V){

  int     k, i, *ai, *jptr ;
  double  *a;
  int found;

  switch (M->T) {
  case SPARSE :
    /* csr_format transpose! 
       donc la matrice arrivant dans cette routine doit
       bel et bien etre la transposee !!! */
    if(M->changed){
      csr_format (&M->S, M->N);
      restore_format (&M->S);
      M->changed = 0 ;
    }
    jptr = (int*) M->S.jptr->array;
    a    = (double*) M->S.a->array;
    ai   = (int*) M->S.ai->array;

    for(i=0; i<M->N; i++){  /* lignes */
      found=0;
      for(k=jptr[i]-1;k<jptr[i+1]-1;k++){ /*colonne */
         if(ai[k]-1==col) { 
	   V[i]=a[k]; found=1; break; 
	 }
	 else if (ai[k]-1 > col) { 
	   break; 
	 }
       }
      if (!found) V[i]=0; 
    printf(" V[%d] = %g \n",i, V[i]); 
    } 
    break;
  case DENSE :
    if(M->notranspose){
	*V = M->F.a[(M->N)*col+row];
    }
    else{
      for(i=0; i<M->N; i++){
	*V = M->F.a[(M->N)*row+col];
      }
    }
    break;
  }
}