Changes between Initial Version and Version 1 of cypress/Programming/Cg


Ignore:
Timestamp:
May 14, 2015 3:17:51 PM (7 years ago)
Author:
cmaggio
Comment:

migrated content from ccs wiki

Legend:

Unmodified
Added
Removed
Modified
  • cypress/Programming/Cg

    v1 v1  
     1= Cg.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5 
     6  solve Ax=b 
     7  Conjugate Gradient Method
     8  General Sparse Matrix A
     9 
     10  solve_cg solve Ax=b with CG method
     11  this returns the number of iterations, or negative if failed.
     12 
     13  Note: A must be symmetric and positive-defined matrix
     14 */
     15#include<stdio.h>
     16 
     17int solve_cg(int dim,          /* dimension */
     18             int nnz,       /* # of non-zeros in the matrix */
     19             double *nzval,    /* array of nonzero values */
     20             int *colidx,      /* array of column indices of the nonzeros */
     21             int *rowptr,      /* the first column of nonzero in each row */
     22             /* rowptr[j] stores the location of nzval[] and colidx[], */
     23             /* which starts row j. This has dim+1 entry that is nnz */
     24             double *x,        /* solition */
     25             double *b,        /* right hand side */
     26             double tol);      /* tolerance */
     27}}}
     28{{{#!C
     29/*
     30  Large Scale Computing
     31 
     32  solve Ax=b 
     33  Conjugate Gradient Method
     34  General Sparse Matrix A
     35 
     36  Note: A must be symmetric and positive-defined matrix
     37 */
     38#include<stdio.h>
     39#include<stdlib.h>
     40#include<math.h>
     41 
     42int solve_cg(int dim,          /* dimension */
     43             int nnz,          /* # of non-zeros in the matrix */
     44             double *nzval,    /* array of nonzero values */
     45             int *colidx,      /* array of column indices of the nonzeros */
     46             int *rowptr,      /* the first column of nonzero in each row */
     47             double *x,        /* solition */
     48             double *b,        /* right hand side */
     49             double tol){      /* tolerance */
     50  double *r;
     51  double *p;
     52  double *Ap;
     53  double rr,pAp,rr1,alpha,beta;
     54  int i,j,count;
     55 
     56  /* Allocate Working Spaces */
     57  r = (double *)calloc(dim, sizeof(double));
     58  p = (double *)calloc(dim, sizeof(double));
     59  Ap = (double *)calloc(dim, sizeof(double));
     60  if ((r == NULL) || (p == NULL) || (Ap == NULL)){
     61    printf("memory allocation failed\n");
     62    return -1;
     63  }
     64 
     65  /* compute r0 */
     66  for (i = 0 ; i < dim ; i++){
     67    r[i] = b[i];
     68    for (j = rowptr[i] ; j < rowptr[i+1] ; j++){
     69      r[i] -= nzval[j] * x[colidx[j]];
     70    }
     71  }
     72 
     73  rr = 0.0;
     74  for (i = 0 ; i < dim ; i++){
     75    p[i] = r[i];    /* p = r */
     76    rr += r[i] * r[i];  /* rr = r.r */
     77  }
     78 
     79  /* cg iteration */
     80  count = 0;
     81  while(rr > tol){   
     82    // Ap = A*p
     83    for (i = 0 ; i < dim ; i++){
     84      Ap[i] = 0.0;
     85      for (j = rowptr[i] ; j < rowptr[i+1] ; j++){
     86        Ap[i] += nzval[j] * p[colidx[j]];
     87      }
     88    }
     89 
     90    // alpha = r.r / p.Ap
     91    pAp = 0.0;
     92    for (i = 0 ; i < dim ; i++){
     93      pAp += p[i] * Ap[i];
     94    }
     95    alpha = rr / pAp;
     96 
     97    //Beta
     98    rr1 = 0.0;
     99    for (i = 0 ; i < dim ; i++){
     100      x[i] += alpha * p[i];
     101      r[i] -= alpha * Ap[i];
     102      rr1 += r[i] * r[i];
     103    } 
     104 
     105    beta = rr1 / rr;
     106    for (i = 0 ; i < dim ; i++){
     107      p[i] = r[i] + beta * p[i];
     108    }
     109 
     110    rr = rr1;
     111    count++;
     112  }
     113 
     114  /* Deallocate Working Spaces */
     115  free(r);
     116  free(p);
     117  free(Ap);
     118 
     119  return count;
     120}
     121}}}