# Diff of /branches/domexper/dudley/src/Util.c

revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC
# Line 97  void Finley_Util_SmallMatSetMult(int len Line 97  void Finley_Util_SmallMatSetMult(int len
97      }      }
98  }  }
99
/* calcultes the LU factorization for a small matrix dimxdim matrix A */
/* TODO: use LAPACK */

int Finley_Util_SmallMatLU(int dim,double* A,double *LU,int* pivot){
double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
int info=0;
/* LAPACK version */
/* memcpy(LU,A,sizeof(douple)); */
/* dgetf2_(&dim,&dim,A,&dim,LU,pivot,&info); */
switch(dim) {
case 1:
D=A[INDEX2(0,0,dim)];
if (ABS(D) >0. ){
LU[INDEX2(0,0,dim)]=1./D;
} else {
info=2;
}
break;

case 2:
A11=A[INDEX2(0,0,dim)];
A12=A[INDEX2(0,1,dim)];
A21=A[INDEX2(1,0,dim)];
A22=A[INDEX2(1,1,dim)];

D = A11*A22-A12*A21;
if (ABS(D) > 0 ){
D=1./D;
LU[INDEX2(0,0,dim)]= A22*D;
LU[INDEX2(1,0,dim)]=-A21*D;
LU[INDEX2(0,1,dim)]=-A12*D;
LU[INDEX2(1,1,dim)]= A11*D;
} else {
info=2;
}
break;

case 3:
A11=A[INDEX2(0,0,dim)];
A21=A[INDEX2(1,0,dim)];
A31=A[INDEX2(2,0,dim)];
A12=A[INDEX2(0,1,dim)];
A22=A[INDEX2(1,1,dim)];
A32=A[INDEX2(2,1,dim)];
A13=A[INDEX2(0,2,dim)];
A23=A[INDEX2(1,2,dim)];
A33=A[INDEX2(2,2,dim)];

D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
if (ABS(D) > 0 ){
D=1./D;
LU[INDEX2(0,0,dim)]=(A22*A33-A23*A32)*D;
LU[INDEX2(1,0,dim)]=(A31*A23-A21*A33)*D;
LU[INDEX2(2,0,dim)]=(A21*A32-A31*A22)*D;
LU[INDEX2(0,1,dim)]=(A13*A32-A12*A33)*D;
LU[INDEX2(1,1,dim)]=(A11*A33-A31*A13)*D;
LU[INDEX2(2,1,dim)]=(A12*A31-A11*A32)*D;
LU[INDEX2(0,2,dim)]=(A12*A23-A13*A22)*D;
LU[INDEX2(1,2,dim)]=(A13*A21-A11*A23)*D;
LU[INDEX2(2,2,dim)]=(A11*A22-A12*A21)*D;
} else {
info=2;
}
break;
default:
info=1;
}
return info;
}

/* solves LUx=b where LU is a LU factorization calculated by an Finley_Util_SmallMatLU call */
void Finley_Util_SmallMatForwardBackwardSolve(int dim ,int nrhs,double* LU,int* pivot,double* x,double* b) {
int i;
switch(dim) {
case 1:
for (i=0;i<nrhs;i++) {
x[INDEX2(0,i,dim)]=LU[0]*b[INDEX2(0,i,dim)];
}
break;
case 2:
for (i=0;i<nrhs;i++) {
x[INDEX2(0,i,dim)]=LU[INDEX2(0,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(0,1,dim)]*b[INDEX2(1,i,dim)];
x[INDEX2(1,i,dim)]=LU[INDEX2(1,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(1,1,dim)]*b[INDEX2(1,i,dim)];
}
break;

case 3:
for (i=0;i<nrhs;i++) {
x[INDEX2(0,i,dim)]=LU[INDEX2(0,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(0,1,dim)]*b[INDEX2(1,i,dim)]+LU[INDEX2(0,2,dim)]*b[INDEX2(2,i,dim)];
x[INDEX2(1,i,dim)]=LU[INDEX2(1,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(1,1,dim)]*b[INDEX2(1,i,dim)]+LU[INDEX2(1,2,dim)]*b[INDEX2(2,i,dim)];
x[INDEX2(2,i,dim)]=LU[INDEX2(2,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(2,1,dim)]*b[INDEX2(1,i,dim)]+LU[INDEX2(2,2,dim)]*b[INDEX2(2,i,dim)];
}
break;
}
return;
}
100  /*    inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */  /*    inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
101  /*    the determinante is returned. */  /*    the determinante is returned. */
102
# Line 204  void Finley_Util_InvertSmallMat(int len, Line 108  void Finley_Util_InvertSmallMat(int len,
108        case 1:        case 1:
109           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
110              D=A[INDEX3(0,0,q,dim,dim)];              D=A[INDEX3(0,0,q,dim,dim)];
111              if (ABS(D) > 0 ){              if (D == 0 ){
det[q]=D;
D=1./D;
invA[INDEX3(0,0,q,dim,dim)]=D;
} else {
112                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
113                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
114                 return;                 return;
115              }              }
116                det[q]=D;
117                D=1./D;
118                invA[INDEX3(0,0,q,dim,dim)]=D;
119           }           }
120           break;           break;
121
# Line 224  void Finley_Util_InvertSmallMat(int len, Line 127  void Finley_Util_InvertSmallMat(int len,
127              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,dim,dim)];
128
129              D = A11*A22-A12*A21;              D = A11*A22-A12*A21;
130              if (ABS(D) > 0 ){              if (D == 0 ){
det[q]=D;
D=1./D;
invA[INDEX3(0,0,q,dim,dim)]= A22*D;
invA[INDEX3(1,0,q,dim,dim)]=-A21*D;
invA[INDEX3(0,1,q,dim,dim)]=-A12*D;
invA[INDEX3(1,1,q,dim,dim)]= A11*D;
} else {
131                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
132                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
133                 return;                 return;
134              }              }
135                det[q]=D;
136                D=1./D;
137                invA[INDEX3(0,0,q,dim,dim)]= A22*D;
138                invA[INDEX3(1,0,q,dim,dim)]=-A21*D;
139                invA[INDEX3(0,1,q,dim,dim)]=-A12*D;
140                invA[INDEX3(1,1,q,dim,dim)]= A11*D;
141           }           }
142           break;           break;
143
# Line 252  void Finley_Util_InvertSmallMat(int len, Line 154  void Finley_Util_InvertSmallMat(int len,
154              A33=A[INDEX3(2,2,q,dim,dim)];              A33=A[INDEX3(2,2,q,dim,dim)];
155
156              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
157              if (ABS(D) > 0 ){              if (D == 0 ){
det[q]    =D;
D=1./D;
invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;
invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;
invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;
invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;
invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;
invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;
invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;
invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;
invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;
} else {
158                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
159                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
160                 return;                 return;
161              }              }
162                det[q]    =D;
163                D=1./D;
164                invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;
165                invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;
166                invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;
167                invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;
168                invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;
169                invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;
170                invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;
171                invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;
172                invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;
173           }           }
174           break;           break;
175
# Line 505  void Finley_copyDouble(int n,double* sou Line 406  void Finley_copyDouble(int n,double* sou
406
407  /*  /*
408   * \$Log\$   * \$Log\$
409   * Revision 1.2  2004/12/14 05:39:31  jgs   * Revision 1.3  2004/12/15 03:48:47  jgs
410   * *** empty log message ***   * *** empty log message ***
411   *   *
* Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross
* a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
*
412   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
413   * initial import of project esys2   * initial import of project esys2
414   *   *

Legend:
 Removed from v.97 changed lines Added in v.100