scilab/modules/sparse/help/en_US/iterativesolvers/qmr.xml

   1 <?xml version="1.0" encoding="UTF-8"?>
   2 <!--
   3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   4  * Copyright (C) XXXX-2008 - INRIA
   5  *
   6  * This file must be used under the terms of the CeCILL.
   7  * This source file is licensed as described in the file COPYING, which
   8  * you should have received as part of this distribution.  The terms
   9  * are also available at
  10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  11  *
  12  -->
  13 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="qmr">
  14     <refnamediv>
  15         <refname>qmr</refname>
  16         <refpurpose>quasi minimal residual method with preconditioning  </refpurpose>
  17     </refnamediv>
  18     <refsynopsisdiv>
  19         <title>Calling Sequence</title>
  20         <synopsis>[x,flag,err,iter,res] = qmr(A,Ap,b,x0,M1,M1p,M2,M2p,maxi,tol)
  21             [x,flag,err,iter,res] = qmr(A,b,x0,M1,M2,maxi,tol)
  22         </synopsis>
  23     </refsynopsisdiv>
  24     <refsection>
  25         <title>Arguments</title>
  26         <variablelist>
  27             <varlistentry>
  28                 <term>A</term>
  29                 <listitem>
  30                     <para>
  31                         matrix of size n-by-n or function.
  32                     </para>
  33                     <itemizedlist>
  34                         <listitem>
  35                             <para>
  36                                 <literal>matrix.</literal>If A is a matrix, it can be
  37                                 dense or sparse
  38                             </para>
  39                         </listitem>
  40                         <listitem>
  41                             <para>
  42                                 <literal>function.</literal>If A is a function which returns <literal>A*x</literal>, it must
  43                                 have the following header :
  44                             </para>
  45                             <programlisting role=""><![CDATA[
  46 function y = A ( x )
  47  ]]></programlisting>
  48                             <para>
  49                                 If A is a function which returns <literal>A*x</literal> or <literal>A'*x</literal> depending t.
  50                                 If <literal>t = "notransp"</literal>, the function returns <literal>A*x</literal>.
  51                                 If <literal>t = "transp"</literal>, the function returns <literal>A'*x</literal>. It must
  52                                 have the following header :
  53                             </para>
  54                             <programlisting role=""><![CDATA[
  55 function y = A ( x, t )
  56  ]]></programlisting>
  57                         </listitem>
  58                     </itemizedlist>
  59                 </listitem>
  60             </varlistentry>
  61             <varlistentry>
  62                 <term>Ap</term>
  63                 <listitem>
  64                     <para>
  65                         function returning <literal>A'*x</literal>. It must have the followinf header :
  66                     </para>
  67                     <programlisting role=""><![CDATA[
  68 function y = Ap ( x )
  69  ]]></programlisting>
  70                 </listitem>
  71             </varlistentry>
  72             <varlistentry>
  73                 <term>b</term>
  74                 <listitem>
  75                     <para>right hand side vector</para>
  76                 </listitem>
  77             </varlistentry>
  78             <varlistentry>
  79                 <term>x0</term>
  80                 <listitem>
  81                     <para>initial guess vector (default: zeros(n,1))</para>
  82                 </listitem>
  83             </varlistentry>
  84             <varlistentry>
  85                 <term>M1</term>
  86                 <listitem>
  87                     <para>
  88                         left preconditioner : matrix or function (In the first case, default: eye(n,n)). If <literal>M1</literal> is a function, she returns either,
  89                     </para>
  90                     <itemizedlist>
  91                         <listitem>
  92                             <para>
  93                                 only <literal>M1*x</literal>
  94                             </para>
  95                         </listitem>
  96                         <para> or
  97                         </para>
  98                         <listitem>
  99                             <para>
 100                                 <literal>M1*x</literal> or <literal>M1'*x</literal> depending <literal>t</literal>.
 101                             </para>
 102                         </listitem>
 103                     </itemizedlist>
 104                 </listitem>
 105             </varlistentry>
 106             <varlistentry>
 107                 <term>M1p</term>
 108                 <listitem>
 109                     <para>
 110                         must only be provided when <literal>M1</literal> is a function returning  <literal>M1*x</literal>.
 111                         In this case <literal>M1p</literal> is the function which returns <literal>M1'*x</literal>.
 112                     </para>
 113                 </listitem>
 114             </varlistentry>
 115             <varlistentry>
 116                 <term>M2</term>
 117                 <listitem>
 118                     <para>
 119                         right preconditioner : matrix or function (In the first case, default: eye(n,n)). If <literal>M2</literal> is a function, she returns either
 120                     </para>
 121                     <itemizedlist>
 122                         <listitem>
 123                             <para>
 124                                 only <literal>M2*x</literal>
 125                             </para>
 126                         </listitem>
 127                         <para> or
 128                         </para>
 129                         <listitem>
 130                             <para>
 131                                 <literal>M2*x</literal> or <literal>M2'*x</literal> depending <literal>t</literal>.
 132                             </para>
 133                         </listitem>
 134                     </itemizedlist>
 135                 </listitem>
 136             </varlistentry>
 137             <varlistentry>
 138                 <term>M2p</term>
 139                 <listitem>
 140                     <para>
 141                         must only be provided when <literal>M2</literal> is a function returning  <literal>M2*x</literal>.
 142                         In this case <literal>M2p</literal> is the function which returns <literal>M2'*x</literal>
 143                     </para>
 144                 </listitem>
 145             </varlistentry>
 146             <varlistentry>
 147                 <term>maxi</term>
 148                 <listitem>
 149                     <para>maximum number of iterations (default: n)
 150                     </para>
 151                 </listitem>
 152             </varlistentry>
 153             <varlistentry>
 154                 <term>tol</term>
 155                 <listitem>
 156                     <para>error tolerance (default: 1000*%eps)</para>
 157                 </listitem>
 158             </varlistentry>
 159             <varlistentry>
 160                 <term>x</term>
 161                 <listitem>
 162                     <para>solution vector</para>
 163                 </listitem>
 164             </varlistentry>
 165             <varlistentry>
 166                 <term>flag</term>
 167                 <listitem>
 168                     <itemizedlist>
 169                         <listitem>
 170                             <para><varname>flag</varname>=0: <code>qmr</code> converged to the desired tolerance within <varname>maxi</varname>
 171                             iterations,</para>
 172                         </listitem>
 173                         <listitem>
 174                             <para><varname>flag</varname>=1: no convergence given <varname>maxi</varname>,</para>
 175                         </listitem>
 176                         <listitem>
 177                             <para><literal>-7 &lt; flag &lt; 0</literal>: A breakdown occurred because one of the scalar quantities calculated during
 178                             <code>qmr</code> was equal to zero.</para>
 179                         </listitem>
 180                     </itemizedlist>
 181                 </listitem>
 182             </varlistentry>
 183             <varlistentry>
 184                 <term>res</term>
 185                 <listitem>
 186                     <para>residual vector</para>
 187                 </listitem>
 188             </varlistentry>
 189             <varlistentry>
 190                 <term>err</term>
 191                 <listitem>
 192                     <para>final residual norm</para>
 193                 </listitem>
 194             </varlistentry>
 195             <varlistentry>
 196                 <term>iter</term>
 197                 <listitem>
 198                     <para>number of iterations performed</para>
 199                 </listitem>
 200             </varlistentry>
 201         </variablelist>
 202     </refsection>
 203     <refsection>
 204         <title>Description</title>
 205         <para>
 206             Solves the linear system <literal>Ax=b</literal> using the Quasi Minimal Residual Method with preconditioning.
 207         </para>
 208     </refsection>
 209     <refsection>
 210         <title>Examples</title>
 211         <programlisting role="example"><![CDATA[
 212         // If A is a matrix
 213 A=[ 94  0   0   0    0   28  0   0   32  0
 214      0   59  13  5    0   0   0   10  0   0
 215      0   13  72  34   2   0   0   0   0   65
 216      0   5   34  114  0   0   0   0   0   55
 217      0   0   2   0    70  0   28  32  12  0
 218      28  0   0   0    0   87  20  0   33  0
 219      0   0   0   0    28  20  71  39  0   0
 220      0   10  0   0    32  0   39  46  8   0
 221      32  0   0   0    12  33  0   8   82  11
 222      0   0   65  55   0   0   0   0   11  100];
 223 b=ones(10,1);
 224 [x,flag,err,iter,res] = qmr(A, b)
 225
 226 [x,flag,err,iter,res] = qmr(A, b, zeros(10,1), eye(10,10), eye(10,10), 10, 1d-12)
 227
 228         // If A is a function
 229         function y = Atimesx(x,t)
 230         A=[ 94  0   0   0    0   28  0   0   32  0
 231      0   59  13  5    0   0   0   10  0   0
 232      0   13  72  34   2   0   0   0   0   65
 233      0   5   34  114  0   0   0   0   0   55
 234      0   0   2   0    70  0   28  32  12  0
 235      28  0   0   0    0   87  20  0   33  0
 236      0   0   0   0    28  20  71  39  0   0
 237      0   10  0   0    32  0   39  46  8   0
 238      32  0   0   0    12  33  0   8   82  11
 239      0   0   65  55   0   0   0   0   11  100];
 240          if (t == 'notransp') then
 241         y = A*x;
 242     elseif (t ==  'transp') then
 243         y = A'*x;
 244     end
 245         endfunction
 246
 247          [x,flag,err,iter,res] = qmr(Atimesx, b)
 248
 249          [x,flag,err,iter,res] = qmr(Atimesx, b, zeros(10,1), eye(10,10), eye(10,10), 10, 1d-12)
 250
 251          // OR
 252
 253          function y = funA(x)
 254         A = [ 94  0   0   0    0   28  0   0   32  0
 255      0   59  13  5    0   0   0   10  0   0
 256      0   13  72  34   2   0   0   0   0   65
 257      0   5   34  114  0   0   0   0   0   55
 258      0   0   2   0    70  0   28  32  12  0
 259      28  0   0   0    0   87  20  0   33  0
 260      0   0   0   0    28  20  71  39  0   0
 261      0   10  0   0    32  0   39  46  8   0
 262      32  0   0   0    12  33  0   8   82  11
 263      0   0   65  55   0   0   0   0   11  100];
 264          y = A*x
 265         endfunction
 266
 267          function y = funAp(x)
 268         A = [ 94  0   0   0    0   28  0   0   32  0
 269      0   59  13  5    0   0   0   10  0   0
 270      0   13  72  34   2   0   0   0   0   65
 271      0   5   34  114  0   0   0   0   0   55
 272      0   0   2   0    70  0   28  32  12  0
 273      28  0   0   0    0   87  20  0   33  0
 274      0   0   0   0    28  20  71  39  0   0
 275      0   10  0   0    32  0   39  46  8   0
 276      32  0   0   0    12  33  0   8   82  11
 277      0   0   65  55   0   0   0   0   11  100];
 278          y = A'*x
 279         endfunction
 280
 281          [x,flag,err,iter,res] = qmr(funA, funAp, b)
 282
 283          [x,flag,err,iter,res] = qmr(funA, funAp, b, zeros(10,1), eye(10,10), eye(10,10), 10, 1d-12)
 284
 285          // If A is a matrix, M1 and M2 are functions
 286          function y = M1timesx(x,t)
 287          M1 = eye(10,10);
 288     if(t=="notransp") then
 289         y = M1*x;
 290     elseif (t=="transp") then
 291         y = M1'*x;
 292     end
 293         endfunction
 294
 295         function y = M2timesx(x,t)
 296          M2 = eye(10,10);
 297     if(t=="notransp") then
 298         y = M2*x;
 299     elseif (t=="transp") then
 300         y = M2'*x;
 301     end
 302         endfunction
 303
 304         [x,flag,err,iter,res] = qmr(A, b, zeros(10,1), M1timesx, M2timesx, 10, 1d-12)
 305
 306         // OR
 307
 308         function y = funM1(x)
 309         M1 = eye(10,10);
 310         y = M1*x;
 311         endfunction
 312
 313         function y = funM1p(x)
 314         M1 = eye(10,10);
 315         y = M1'*x;
 316         endfunction
 317
 318         function y = funM2(x)
 319         M2 = eye(10,10);
 320         y = M2*x;
 321         endfunction
 322
 323         function y = funM2p(x)
 324         M2 = eye(10,10);
 325         y = M2'*x;
 326         endfunction
 327
 328         [x,flag,err,iter,res] = qmr(A, b, zeros(10,1), funM1, funM1p, funM2, funM2p, 10, 1d-12)
 329
 330         // If A, M1, M2 are functions
 331         [x,flag,err,iter,res] = qmr(funA, funAp, b, zeros(10,1), funM1, funM1p, funM2, funM2p, 10, 1d-12)
 332         [x,flag,err,iter,res] = qmr(Atimesx, b, zeros(10,1), M1timesx, M2timesx, 10, 1d-12)
 333
 334  ]]></programlisting>
 335     </refsection>
 336     <refsection role="see also">
 337         <title>See Also</title>
 338         <simplelist type="inline">
 339             <member>
 340                 <link linkend="gmres">gmres</link>
 341             </member>
 342             <member>
 343                 <link linkend="pcg">pcg</link>
 344             </member>
 345         </simplelist>
 346     </refsection>
 347     <refsection>
 348         <title>History</title>
 349         <revhistory>
 350             <revision>
 351                 <revnumber>5.4.0</revnumber>
 352                 <revdescription>
 353                     Calling qmr(A, Ap) is deprecated. qmr(A) should be used instead. The following function is an example :
 354                     <programlisting role=""><![CDATA[
 355 function y = A ( x, t )
 356 Amat = eye(2,2);
 357 if ( t== "notransp") then
 358 y = Amat*x;
 359 elseif (t == "transp") then
 360 y = Amat'*x;
 361 end
 362 endfunction
 363  ]]></programlisting>
 364                 </revdescription>
 365             </revision>
 366         </revhistory>
 367     </refsection>
 368 </refentry>