← Back to branch summary

~ubuntu-branches/ubuntu/warty/petsc/warty

« back to all changes in this revision

Viewing changes to src/sles/ksp/impls/cr/cr.c.html

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2004-06-07 13:41:43 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040607134143-92p586zrauvie0le
Tags: 2.2.0-2
http://bugs.debian.org/249617
* Upstream patch level 2.
* New PETSC_BOPT_EXTRA option for different BOPT and lib names, with _c++
  symlinks only for plain and single (closes: #249617).
* New DEBIAN_DIST=contrib option to link with hypre, parmetis (closes:
  #249619).
* Combined petsc-c and petsc-fortran substvars into petsc-compilers.
* Extra quote in -dev prerm eliminates "too many arguments" problem.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/sles/ksp/impls/cr/cr.c.html
1
 
<center><a href="cr.c">Actual source code: cr.c</a></center><br>
2
 
 
3
 
<html>
4
 
<head>
5
 
<title></title>
6
 
<meta name="generator" content="c2html 0.9.1">
7
 
<meta name="date" content="2002-05-31T16:26:13+00:00">
8
 
</head>
9
 
 
10
 
<body bgcolor="#FFFFFF">
11
 
<pre width="80"><a name="line1">  1: </a><font color="#B22222">/*$Id: cr.c,v 1.64 2001/08/07 03:03:49 balay Exp $*/</font>
12
 
 
13
 
<a name="line3">  3: </a><font color="#B22222">/*                       </font>
14
 
<a name="line4">  4: </a><font color="#B22222">           This implements Preconditioned Conjugate Residuals.       </font>
15
 
<a name="line5">  5: </a><font color="#B22222">*/</font>
16
 
<a name="line6"> 6: </a> #include <A href="../../../../../src/sles/ksp/kspimpl.h.html">src/sles/ksp/kspimpl.h</A>
17
 
 
18
 
<a name="line8">  8: </a><strong><font color="#4169E1"><a name="KSPSetUp_CR"></a>static int KSPSetUp_CR(<A href="../../../../../docs/manualpages/KSP/KSP.html#KSP">KSP</A> ksp)</font></strong>
19
 
<a name="line9">  9: </a>{
20
 
 
21
 
<a name="line13"> 13: </a>  <font color="#4169E1">if</font> (ksp-&gt;pc_side == PC_RIGHT) {<A href="../../../../../docs/manualpages/Sys/SETERRQ.html#SETERRQ">SETERRQ</A>(2,<font color="#666666">"no right preconditioning for KSPCR"</font>);}
22
 
<a name="line14"> 14: </a>  <font color="#4169E1">else</font> <font color="#4169E1">if</font> (ksp-&gt;pc_side == PC_SYMMETRIC) {<A href="../../../../../docs/manualpages/Sys/SETERRQ.html#SETERRQ">SETERRQ</A>(2,<font color="#666666">"no symmetric preconditioning for KSPCR"</font>);}
23
 
<a name="line15"> 15: </a>  KSPDefaultGetWork(ksp,6);
24
 
<a name="line16"> 16: </a>  <font color="#4169E1">return</font>(0);
25
 
<a name="line17"> 17: </a>}
26
 
 
27
 
<a name="line19"> 19: </a><strong><font color="#4169E1"><a name="KSPSolve_CR"></a>static int  KSPSolve_CR(<A href="../../../../../docs/manualpages/KSP/KSP.html#KSP">KSP</A> ksp,int *its)</font></strong>
28
 
<a name="line20"> 20: </a>{
29
 
<a name="line21"> 21: </a>  int          i = 0, maxit, ierr;
30
 
<a name="line22"> 22: </a>  <A href="../../../../../docs/manualpages/Mat/MatStructure.html#MatStructure">MatStructure</A> pflag;
31
 
<a name="line23"> 23: </a>  PetscReal    dp;
32
 
<a name="line24"> 24: </a>  <A href="../../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  ai, bi;
33
 
<a name="line25"> 25: </a>  <A href="../../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  apq,btop, bbot, mone = -1.0;
34
 
<a name="line26"> 26: </a>  <A href="../../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>          X,B,R,RT,P,AP,ART,Q;
35
 
<a name="line27"> 27: </a>  <A href="../../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A>          Amat, Pmat;
36
 
 
37
 
 
38
 
<a name="line31"> 31: </a>  maxit   = ksp-&gt;max_it;
39
 
<a name="line32"> 32: </a>  X       = ksp-&gt;vec_sol;
40
 
<a name="line33"> 33: </a>  B       = ksp-&gt;vec_rhs;
41
 
<a name="line34"> 34: </a>  R       = ksp-&gt;work[0];
42
 
<a name="line35"> 35: </a>  RT      = ksp-&gt;work[1];
43
 
<a name="line36"> 36: </a>  P       = ksp-&gt;work[2];
44
 
<a name="line37"> 37: </a>  AP      = ksp-&gt;work[3];
45
 
<a name="line38"> 38: </a>  ART     = ksp-&gt;work[4];
46
 
<a name="line39"> 39: </a>  Q       = ksp-&gt;work[5];
47
 
 
48
 
<a name="line41"> 41: </a>  <font color="#B22222">/* R is the true residual norm, RT is the preconditioned residual norm */</font>
49
 
<a name="line42"> 42: </a>  <A href="../../../../../docs/manualpages/PC/PCGetOperators.html#PCGetOperators">PCGetOperators</A>(ksp-&gt;B,&amp;Amat,&amp;Pmat,&amp;pflag);
50
 
<a name="line43"> 43: </a>  <font color="#4169E1">if</font> (!ksp-&gt;guess_zero) {
51
 
<a name="line44"> 44: </a>    KSP_MatMult(ksp,Amat,X,R);     <font color="#B22222">/*   R &lt;- A*X           */</font>
52
 
<a name="line45"> 45: </a>    <A href="../../../../../docs/manualpages/Vec/VecAYPX.html#VecAYPX">VecAYPX</A>(&amp;mone,B,R);            <font color="#B22222">/*   R &lt;- B-R == B-A*X  */</font>
53
 
<a name="line46"> 46: </a>  } <font color="#4169E1">else</font> {
54
 
<a name="line47"> 47: </a>    <A href="../../../../../docs/manualpages/Vec/VecCopy.html#VecCopy">VecCopy</A>(B,R);                  <font color="#B22222">/*   R &lt;- B (X is 0)    */</font>
55
 
<a name="line48"> 48: </a>  }
56
 
<a name="line49"> 49: </a>  KSP_PCApply(ksp,ksp-&gt;B,R,P);     <font color="#B22222">/*   P   &lt;- B*R         */</font>
57
 
<a name="line50"> 50: </a>  KSP_MatMult(ksp,Amat,P,AP);      <font color="#B22222">/*   AP  &lt;- A*P         */</font>
58
 
<a name="line51"> 51: </a>  <A href="../../../../../docs/manualpages/Vec/VecCopy.html#VecCopy">VecCopy</A>(P,RT);                   <font color="#B22222">/*   RT  &lt;- P           */</font>
59
 
<a name="line52"> 52: </a>  <A href="../../../../../docs/manualpages/Vec/VecCopy.html#VecCopy">VecCopy</A>(AP,ART);                 <font color="#B22222">/*   ART &lt;- AP          */</font>
60
 
<a name="line53"> 53: </a>  ierr   = <A href="../../../../../docs/manualpages/Vec/VecDot.html#VecDot">VecDot</A>(RT,ART,&amp;btop);          <font color="#B22222">/*   (RT,ART)           */</font>
61
 
<a name="line54"> 54: </a>  <font color="#4169E1">if</font> (ksp-&gt;normtype == KSP_PRECONDITIONED_NORM) {
62
 
<a name="line55"> 55: </a>    <A href="../../../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(RT,NORM_2,&amp;dp);        <font color="#B22222">/*   dp &lt;- RT'*RT       */</font>
63
 
<a name="line56"> 56: </a>  } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (ksp-&gt;normtype == KSP_UNPRECONDITIONED_NORM) {
64
 
<a name="line57"> 57: </a>    <A href="../../../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(R,NORM_2,&amp;dp);         <font color="#B22222">/*   dp &lt;- R'*R         */</font>
65
 
<a name="line58"> 58: </a>  } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (ksp-&gt;normtype == KSP_NATURAL_NORM) {
66
 
<a name="line59"> 59: </a>    dp = sqrt(PetscAbsScalar(btop));                    <font color="#B22222">/* dp = sqrt(R,AR)      */</font>
67
 
<a name="line60"> 60: </a>  }
68
 
<a name="line61"> 61: </a>  (*ksp-&gt;converged)(ksp,0,dp,&amp;ksp-&gt;reason,ksp-&gt;cnvP);
69
 
<a name="line62"> 62: </a>  <font color="#4169E1">if</font> (ksp-&gt;reason) {*its = 0; <font color="#4169E1">return</font>(0);}
70
 
<a name="line63"> 63: </a>  KSPMonitor(ksp,0,dp);
71
 
<a name="line64"> 64: </a>  PetscObjectTakeAccess(ksp);
72
 
<a name="line65"> 65: </a>  ksp-&gt;its = 0;
73
 
<a name="line66"> 66: </a>  ksp-&gt;rnorm              = dp;
74
 
<a name="line67"> 67: </a>  PetscObjectGrantAccess(ksp);
75
 
<a name="line68"> 68: </a>  KSPLogResidualHistory(ksp,dp);
76
 
 
77
 
<a name="line70"> 70: </a>  <font color="#4169E1">for</font> ( i=0; i&lt;maxit; i++) {
78
 
<a name="line71"> 71: </a>    ierr   = KSP_PCApply(ksp,ksp-&gt;B,AP,Q);<font color="#B22222">/*   Q &lt;- B* AP          */</font>
79
 
<a name="line72"> 72: </a>                                                        <font color="#B22222">/* Step 3                */</font>
80
 
 
81
 
<a name="line74"> 74: </a>    ierr   = <A href="../../../../../docs/manualpages/Vec/VecDot.html#VecDot">VecDot</A>(AP,Q,&amp;apq);
82
 
<a name="line75"> 75: </a>    ai = btop/apq;                                      <font color="#B22222">/* ai = (RT,ART)/(AP,Q)  */</font>
83
 
 
84
 
<a name="line77"> 77: </a>    ierr   = <A href="../../../../../docs/manualpages/Vec/VecAXPY.html#VecAXPY">VecAXPY</A>(&amp;ai,P,X);            <font color="#B22222">/*   X   &lt;- X + ai*P     */</font>
85
 
<a name="line78"> 78: </a>    ai     = -ai;
86
 
<a name="line79"> 79: </a>    ierr   = <A href="../../../../../docs/manualpages/Vec/VecAXPY.html#VecAXPY">VecAXPY</A>(&amp;ai,Q,RT);           <font color="#B22222">/*   RT  &lt;- RT - ai*Q    */</font>
87
 
<a name="line80"> 80: </a>    ierr   = KSP_MatMult(ksp,Amat,RT,ART);<font color="#B22222">/*   ART &lt;-   A*RT       */</font>
88
 
<a name="line81"> 81: </a>    bbot = btop;
89
 
<a name="line82"> 82: </a>    ierr   = <A href="../../../../../docs/manualpages/Vec/VecDot.html#VecDot">VecDot</A>(RT,ART,&amp;btop);
90
 
 
91
 
<a name="line84"> 84: </a>    <font color="#4169E1">if</font> (ksp-&gt;normtype == KSP_PRECONDITIONED_NORM) {
92
 
<a name="line85"> 85: </a>      <A href="../../../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(RT,NORM_2,&amp;dp);      <font color="#B22222">/*   dp &lt;- || RT ||      */</font>
93
 
<a name="line86"> 86: </a>    } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (ksp-&gt;normtype == KSP_NATURAL_NORM) {
94
 
<a name="line87"> 87: </a>      dp = sqrt(PetscAbsScalar(btop));                  <font color="#B22222">/* dp = sqrt(R,AR)       */</font>
95
 
<a name="line88"> 88: </a>    } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (ksp-&gt;normtype == KSP_NO_NORM) {
96
 
<a name="line89"> 89: </a>      dp = 0.0;
97
 
<a name="line90"> 90: </a>    } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (ksp-&gt;normtype == KSP_UNPRECONDITIONED_NORM) {
98
 
<a name="line91"> 91: </a>      <A href="../../../../../docs/manualpages/Vec/VecAXPY.html#VecAXPY">VecAXPY</A>(&amp;ai,AP,R);           <font color="#B22222">/*   R   &lt;- R - ai*AP    */</font>
99
 
<a name="line92"> 92: </a>      <A href="../../../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(R,NORM_2,&amp;dp);       <font color="#B22222">/*   dp &lt;- R'*R          */</font>
100
 
<a name="line93"> 93: </a>    } <font color="#4169E1">else</font> {
101
 
<a name="line94"> 94: </a>      <A href="../../../../../docs/manualpages/Sys/SETERRQ1.html#SETERRQ1">SETERRQ1</A>(1,<font color="#666666">"<A href="../../../../../docs/manualpages/KSP/KSPNormType.html#KSPNormType">KSPNormType</A> of %d not supported"</font>,(int)ksp-&gt;normtype);
102
 
<a name="line95"> 95: </a>    }
103
 
 
104
 
<a name="line97"> 97: </a>    PetscObjectTakeAccess(ksp);
105
 
<a name="line98"> 98: </a>    ksp-&gt;its++;
106
 
<a name="line99"> 99: </a>    ksp-&gt;rnorm = dp;
107
 
<a name="line100">100: </a>    PetscObjectGrantAccess(ksp);
108
 
 
109
 
<a name="line102">102: </a>    KSPLogResidualHistory(ksp,dp);
110
 
<a name="line103">103: </a>    KSPMonitor(ksp,i+1,dp);
111
 
<a name="line104">104: </a>    (*ksp-&gt;converged)(ksp,i+1,dp,&amp;ksp-&gt;reason,ksp-&gt;cnvP);
112
 
<a name="line105">105: </a>    <font color="#4169E1">if</font> (ksp-&gt;reason) <font color="#4169E1">break</font>;
113
 
 
114
 
<a name="line107">107: </a>    bi = btop/bbot;
115
 
<a name="line108">108: </a>    <A href="../../../../../docs/manualpages/Vec/VecAYPX.html#VecAYPX">VecAYPX</A>(&amp;bi,RT,P);              <font color="#B22222">/*   P &lt;- RT + Bi P     */</font>
116
 
<a name="line109">109: </a>    <A href="../../../../../docs/manualpages/Vec/VecAYPX.html#VecAYPX">VecAYPX</A>(&amp;bi,ART,AP);            <font color="#B22222">/*   AP &lt;- ART + Bi AP  */</font>
117
 
<a name="line110">110: </a>  }
118
 
<a name="line111">111: </a>  <font color="#4169E1">if</font> (i == maxit) {
119
 
<a name="line112">112: </a>    ksp-&gt;reason =  KSP_DIVERGED_ITS;
120
 
<a name="line113">113: </a>    i--;
121
 
<a name="line114">114: </a>  }
122
 
<a name="line115">115: </a>  *its = i + 1;
123
 
<a name="line116">116: </a>  <font color="#4169E1">return</font>(0);
124
 
<a name="line117">117: </a>}
125
 
 
126
 
 
127
 
<a name="line120">120: </a><strong><font color="#4169E1"><a name="KSPCreate_CR"></a>EXTERN_C_BEGIN</font></strong>
128
 
<a name="line121">121: </a><strong><font color="#4169E1">int KSPCreate_CR(<A href="../../../../../docs/manualpages/KSP/KSP.html#KSP">KSP</A> ksp)</font></strong>
129
 
<a name="line122">122: </a>{
130
 
<a name="line124">124: </a>  ksp-&gt;pc_side                   = PC_LEFT;
131
 
<a name="line125">125: </a>  ksp-&gt;ops-&gt;setup                = KSPSetUp_CR;
132
 
<a name="line126">126: </a>  ksp-&gt;ops-&gt;solve                = KSPSolve_CR;
133
 
<a name="line127">127: </a>  ksp-&gt;ops-&gt;destroy              = KSPDefaultDestroy;
134
 
<a name="line128">128: </a>  ksp-&gt;ops-&gt;buildsolution        = KSPDefaultBuildSolution;
135
 
<a name="line129">129: </a>  ksp-&gt;ops-&gt;buildresidual        = KSPDefaultBuildResidual;
136
 
<a name="line130">130: </a>  ksp-&gt;ops-&gt;setfromoptions       = 0;
137
 
<a name="line131">131: </a>  ksp-&gt;ops-&gt;view                 = 0;
138
 
<a name="line132">132: </a>  <font color="#4169E1">return</font>(0);
139
 
<a name="line133">133: </a>}
140
 
<a name="line134">134: </a>EXTERN_C_END
141
 
</pre>
142
 
</body>
143
 
 
144
 
</html>