1 /* 2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT, 3 * Applied Mathematics, Norway. 4 * 5 * Contact information: E-mail: tor.dokken@sintef.no 6 * SINTEF ICT, Department of Applied Mathematics, 7 * P.O. Box 124 Blindern, 8 * 0314 Oslo, Norway. 9 * 10 * This file is part of SISL. 11 * 12 * SISL is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU Affero General Public License as 14 * published by the Free Software Foundation, either version 3 of the 15 * License, or (at your option) any later version. 16 * 17 * SISL is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 * GNU Affero General Public License for more details. 21 * 22 * You should have received a copy of the GNU Affero General Public 23 * License along with SISL. If not, see 24 * <http://www.gnu.org/licenses/>. 25 * 26 * In accordance with Section 7(b) of the GNU Affero General Public 27 * License, a covered work must retain the producer line in every data 28 * file that is created or manipulated using SISL. 29 * 30 * Other Usage 31 * You can be released from the requirements of the license by purchasing 32 * a commercial license. Buying such a license is mandatory as soon as you 33 * develop commercial activities involving the SISL library without 34 * disclosing the source code of your own applications. 35 * 36 * This file may be used in accordance with the terms contained in a 37 * written agreement between you and SINTEF ICT. 38 */ 39 40 #include "sisl-copyright.h" 41 42 /* 43 * 44 * $Id: sh1924.c,v 1.2 2001-03-19 15:59:06 afr Exp $ 45 * 46 */ 47 48 #define SH1924 49 50 #include "sislP.h" 51 52 #if defined(SISLNEEDPROTOTYPES) 53 void sh1924(double * ea,double * eb,int in,int ik,int idim,int * nstart,int * jstat)54 sh1924(double *ea,double *eb,int in,int ik,int idim, 55 int *nstart,int *jstat) 56 #else 57 void sh1924(ea,eb,in,ik,idim,nstart,jstat) 58 double *ea; 59 double *eb; 60 int in; 61 int ik; 62 int idim; 63 int *nstart; 64 int *jstat; 65 #endif 66 /* 67 ********************************************************************* 68 * 69 * PURPOSE : To solve idim linear systems of equations A*X=eb, given 70 * the right-hand-side eb and the Cholesky factorization of 71 * the matrix A (i.e. A is assumed to be symmetric and 72 * positive definite). It is also assumed that A has at most 73 * 2*ik-1 nonzero elements in each row, and since this 74 * structure is not destroyed by the Cholesky factorization, 75 * it is sufficient to store a maximum of ik elements per 76 * row, see the subroutine sh1923. The right-hand-side eb has 77 * idim components each of length in. The output of the 78 * routine is the solution fo the linear system which 79 * overwrites eb. 80 * 81 * 82 * INPUT : ea - Real array of dimension (in*ik) containing the 83 * nonzero elements of the Cholesky factorization. 84 * eb - Real array of dimension (in*idim) containing 85 * the right-hand-side(s) of the linear systems. 86 * in - The dimension of the linear systems, i.e. the number 87 * of rows in ea. 88 * ik - The maximum number of different nonzero elements 89 * in each row of the Cholesky factorization and 90 * therefore the number of columns of ea. 91 * idim - The number of different right-hand-sides. 92 * nstart - Integer array of dimension (in) containing pointers 93 * to the first nonzero element of each row of ea. 94 * 95 * 96 * OUTPUT : eb - The solution of the idim linear systems. 97 * jstat - status messages 98 * > 0 : warning 99 * = 0 : ok 100 * < 0 : error 101 * 102 * 103 * METHOD : The linear system is solved in the usual way, first a 104 * forward substition and the a back substitution. 105 * 106 * 107 * REFERENCES : Any book on general numerical analysis or numerical 108 * linear algebra. 109 * 110 * 111 * USE : 112 * 113 *- 114 * CALLS : 115 * 116 * WRITTEN BY : Vibeke Skytt, SI, 05.92, on the basis of a routine 117 * written by Tom Lyche and Knut Moerken, 12.85. 118 * 119 ********************************************************************* 120 */ 121 { 122 int ki,kj,kr; /* Counters. */ 123 int ki1,kih,kim; /* Counters. */ 124 int kjs,kjh; /* Pointers into matrix. */ 125 int kik1 = ik-1; /* Order minus one. */ 126 double thelp; /* Help variable. */ 127 double *ssum=SISL_NULL; /* Help array. */ 128 129 /* Allocate scratch for help array. */ 130 131 if ((ssum = new0array(idim,DOUBLE)) == SISL_NULL) goto err101; 132 133 /* Forward substitution. */ 134 135 for (ki=0; ki<in; ki++) 136 { 137 ki1 = ki - 1; 138 for (kr=0; kr<idim; kr++) ssum[kr] = (double)0.0; 139 140 /* kjs points to the firs nonzero element of row ki of ea, and 141 kjh points to the corresponding element in the underlying 142 matrix. */ 143 144 for (kjs=nstart[ki], kjh=ki+kjs-ik+1, kj=kjs; 145 kj<kik1; kjh++, kj++) 146 { 147 thelp = ea[ki*ik+kj]; 148 for (kr=0; kr<idim; kr++) 149 ssum[kr] += thelp*eb[kjh*idim+kr]; 150 } 151 152 /* Check if the linear system is singular. */ 153 154 if (DEQUAL(ea[ki*ik+kik1],DZERO)) goto err106; 155 156 thelp = (double)1.0/ea[ki*ik+kik1]; 157 for (kr=0; kr<idim; kr++) 158 eb[ki*idim+kr] = (eb[ki*idim+kr] - ssum[kr])*thelp; 159 } 160 161 /* Back substitution. */ 162 163 for (ki=in-1, ki1=in-1, kih=0; 164 kih<in; ki--, kih++) 165 { 166 /* Compute the index of the last nonzero element in row ki of the 167 transpose of the Cholesky factorization. The integer ik-ki1+ki 168 gives the position of element no (ki1*ik+ki) in row ki1 of ea. 169 ki1 is reduced until the first nonzero element or row ki1 is 170 ik-ki1+ki. */ 171 172 for (;;ki1--) 173 if (nstart[ki1] < ik-ki1+ki) break; 174 175 for (kr=0; kr<idim; kr++) ssum[kr] = (double)0.0; 176 177 /* Calculate eb[.*ik+kik1]. */ 178 179 for (kjs=ki+1, kim=ik-kjs+ki-1, kj=kjs; 180 kj<=ki1; kim--, kj++) 181 { 182 thelp = ea[kj*ik+kim]; 183 for (kr=0; kr<idim; kr++) 184 ssum[kr] += thelp*eb[kj*idim+kr]; 185 } 186 187 thelp = (double)1.0/ea[ki*ik+ik-1]; 188 for (kr=0; kr<idim; kr++) 189 eb[ki*idim+kr] = (eb[ki*idim+kr] - ssum[kr])*thelp; 190 } 191 192 193 /* The linear system is solved. */ 194 195 *jstat = 0; 196 goto out; 197 198 /* Error in space allocation. */ 199 200 err101: *jstat = -101; 201 goto out; 202 203 /* Singular matrix. */ 204 205 err106: *jstat = -106; 206 goto out; 207 208 out: 209 /* Free scratch used for local array. */ 210 211 if (ssum != SISL_NULL) freearray(ssum); 212 213 return; 214 } 215 216