Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
givens.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #include "fasp.h"
9 #include "fasp_functs.h"
10 
11 /*---------------------------------*/
12 /*-- Public Functions --*/
13 /*---------------------------------*/
14 
28 void fasp_aux_givens (const REAL beta,
29  dCSRmat *H,
30  dvector *y,
31  REAL *tmp)
32 {
33  const INT Hsize=H->row;
34  INT i, j, istart, idiag, ip1start;
35  REAL h0,h1,r,c,s,tempi,tempip1,sum;
36 
37  tmp[0]=beta;
38  memset(&tmp[1], 0x0, sizeof(REAL)*(Hsize-1));
39 
40  for (i=0;i<Hsize-1;++i) {
41  istart=H->IA[i];
42  ip1start=H->IA[i+1];
43  if (i==0) idiag=istart;
44  else idiag=istart+1;
45 
46  h0=H->val[idiag]; // h0=H[i][i]
47  h1=H->val[H->IA[i+1]]; // h1=H[i+1][i]
48  r=sqrt(h0*h0+h1*h1);
49  c=h0/r; s=h1/r;
50 
51  for (j=idiag;j<ip1start;++j) {
52  tempi=H->val[j];
53  tempip1=H->val[ip1start+(j-idiag)];
54  H->val[j]=c*tempi+s*tempip1;
55  H->val[ip1start+(j-idiag)]=c*tempip1-s*tempi;
56  }
57 
58  tempi=c*tmp[i]+s*tmp[i+1];
59  tempip1=c*tmp[i+1]-s*tmp[i];
60 
61  tmp[i]=tempi; tmp[i+1]=tempip1;
62  }
63 
64  for (i=Hsize-2;i>=0;--i) {
65  sum=tmp[i];
66  istart=H->IA[i];
67  if (i==0) idiag=istart;
68  else idiag=istart+1;
69 
70  for (j=Hsize-2;j>i;--j) sum-=H->val[idiag+j-i]*y->val[j];
71 
72  y->val[i]=sum/H->val[idiag];
73  }
74 
75 }
76 
77 /*---------------------------------*/
78 /*-- End of File --*/
79 /*---------------------------------*/
#define REAL
Definition: fasp.h:67
REAL * val
nonzero entries of A
Definition: fasp.h:166
REAL * val
actual vector entries
Definition: fasp.h:348
Vector with n entries of REAL type.
Definition: fasp.h:342
#define INT
Definition: fasp.h:64
void fasp_aux_givens(const REAL beta, dCSRmat *H, dvector *y, REAL *tmp)
Perform Givens rotations to compute y |beta*e_1- H*y|.
Definition: givens.c:28
Main header file for FASP.
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160